Development and Homeostasis of the Nervous System - SC18139
pachnisv Lab - Tiffany Heanue
true
Compiled: November 18, 2020
Abstract
" Glial cells of the vertebrate enteric nervous system ENS have diverse proposed functional properties, ranging from neural progenitors to support cells for neurons F our subtypes of ENS glial have been identified in mouse by morphological and spatial criter i a, however molecular identifiers have remained elusive A variety of projects are underway in the lab to molecularly characterize mouse glial cells, including bulk and single cell transcriptomics As part of parallel projects examining the development and function of the ENS in zebrafish, a project to characterize zebrafish ENS glial cells was undertaken Surprisingly, none of the established ENS glial cell markers are expressed in the developing or mature zebrafish ENS However, while the majority of the ENS lineage marked by the sox 10 cre nlm-Cherry fl transgenic are neurons, a proportion of these are non neuronal cells These non-neuronal cells are labelled by the her4 1GFP transgenic, a tool which also labels radial glial cells of the zebrafish central nervous system These findings suggest that this non-neuronal population may be glial cell-like T he orthologue of her4 in mouse is Hes5 Intriguingly, using a Hes5GFP transgenic, we have now shown that Hes5GFP cells comprise a subpopulation of mammalian ENS glial cells, constituting a first molecular hallmark of a mammalian ENS glial subtype This underscores the power of our parallel studies in different species We would like to characterize the zebrafish ENS more extensively, including this glial-like population Although we are unable to effectively dissociate the zebrafish gut smooth muscle layer s to isolate the her4 1GFP cells embedded within, we can readily isolate nuclei of the ENS lineage by virtue of the nuclear Cherry of the sox 10 cre nlm-Cherry fl transgenic This resource provides a tool to generate a transcriptomic profile of the zebrafish ENS S ingle nuclei transcriptomics will enable us to establish clear profiles of distinct neuronal and glial-like populations Comparisons of transcriptomic data from zebrafish with analogous mouse data generated by our lab will benefit the future studies in both species Experimental Approach Please include details such as the number of biological and or technical replicates Nuclei from the enteric nervous system ENS lineage can be isolated from adult zebrafish gut tissues by virtue of lineage-specific nuclear fluorescent reporter expression sox 10 cre nlm-Cherry fl and FACS isolation as performed for our bulk analysis We plan to isolate 10, 000 nuclei to load the 10x QC of nuclei can take advantage of the fact that isolated nuclei are DAPI and cherry, and thus visible by fluorescence microscopy This experiment is a follow up on our previous bulk nuclei isolation The single nuclei strategy will enable characterization of the heterogeneous lineage labelled by sox 10 cre nlm-Cherry fl and delineation of the distinct p rofiles of neurons and glial-like cells, and subtypes therein We hope to compare mRNA expression analysis across bulk nuclei and single nuclei mRNAseq data sets, but also to compare to equivalent populations that our lab has profiled in the mouse bulk and single cell mRNAseq full genotype Tg-4725sox 10 cre _ba74 ßactin loxp-stop-loxp-nlmCherry Who wil l perform Bioinformatics Analyses If you are not using a member of the core Bioinformatics team at Crick, please provide some details of how you will perform the analyses We hope to have bioinformatics analysis performed by the relevant member s of the core Bioinformatics team We would like to perform bioinformatics analysis across our zebrafish data sets bulk and single mRNAseq, but also to compare to equivalent populations that our lab has profiled in the mouse both bulk and single cell mRNAseq Indicative Costing ASF Use Parameter Value Instrument Chromium-10X Experiment type 10X-3-RNA-Seq Experiment kit 10X_3prime_mRNA Total number of cells to capture 4000 Number of cells to input given capture rate 5600 Number of Channels 1 Sequencer HiSeq4000 Read type PE100 Reads per cell 100000 Total Reads 400000000 Number of HiSeq lanes MiSeq flowcells 1 Chromium-10X cost 1100 00 HiSeq4000 cost 1500 00 Staff and overheads cost 520 00 Total cost 3120 00"
Data Preparation
Set Environment and load R
## Get interactive session ##
# srun --time=08:00:00 --mem=40G -p int --pty bash
# module purge;source /camp/stp/babs/working/software/modulepath_new_software_tree_2018-08-13;module load pandoc/2.2.3.2-foss-2016b;ml R/3.6.0-foss-2016b-BABS;R;
# sbatch --time=08:00:00 --wrap "module purge;source /camp/stp/babs/working/software/modulepath_new_software_tree_2018-08-13;module load pandoc/2.2.3.2-foss-2016b;ml R/3.6.0-foss-2016b-BABS;Rscript runB.r" --job-name="rB" --mem=42G -o rB.slurm >> commands.txt
# sbatch --time=12:00:00 --wrap "module purge;source /camp/stp/babs/working/software/modulepath_new_software_tree_2018-08-13;module load pandoc/2.2.3.2-foss-2016b;ml R/3.6.0-foss-2016b-BABS;Rscript runB.r" --job-name="rB" -p hmem --mem=300G -o rB.slurm >> commands.txt
Set Parameters
library(dplyr)
library(Seurat)
library(ggplot2)
library(tidyverse)
library(tidyr)
library(knitr)
## Setup plot collection object
VersionPdfExt <- paste0(".V", gsub("-", "", Sys.Date()), ".pdf")
if (dir.exists("/Volumes/babs/working/boeings/")){
hpc.mount <- "/Volumes/babs/working/boeings/"
} else if (dir.exists("Y:/working/boeings/")){
hpc.mount <- "Y:/working/boeings/"
} else if (dir.exists("/camp/stp/babs/working/boeings/")){
hpc.mount <- "/camp/stp/babs/working/boeings/"
} else {
hpc.mount <- ""
}
# source(
# paste0(
# hpc.mount,
# "Stefan/protocol_files/github/boeings/packages/packageSourceCode/SBwebtools.pckg.r"
# )
# )
# source(
# paste0(
# hpc.mount,
# "Stefan/protocol_files/github/boeings/packages/scTools/scTools.r"
# )
# )
source("assets/scTools.r")
source("assets/SBwebtools.pckg.r")
####
FN <- paste0(hpc.mount, "Projects/reference_data/documentation/BC.parameters.txt")
dbTable <- read.delim(
FN,
sep = "\t",
stringsAsFactors = F
)
db.pwd <- as.vector(dbTable[1,1])
if (length(.libPaths()) > 2){
.libPaths(.libPaths()[2:3])
}
ObioFN <- paste0("../", list.files("..")[grep(".bioLOGIC.Robj", list.files(".."))])
load(ObioFN)
# checkFile = paste0(
# Obio@parameterList$project_id,
# ".bioLOGIC.Robj"
# )
#
# if (ObioFN != checkFile){
# exit()
# }
Obio <- setMountingPoint(Obio)
Obio <- setAnalysisPaths(Obio)
Obio <- setCrickGenomeAndGeneNameTable(Obio)
Obio <- createAnalysisFolders(
Obio,
baseDir="/camp/stp/babs/working/boeings/Projects/",
localBaseDir = paste0(hpc.mount, "Projects/")
)
Obio <- setDataBaseParameters(Obio)
Obio <- addGeneAnnotation(Obio)
Obio@parameterList[["reportFigDir"]] <- paste0(Obio@parameterList$localWorkDir,Obio@parameterList$project_id, "/report_figures/")
## Create outputfolders ##
if (!dir.exists(paste0(Obio@parameterList$localWorkDir,Obio@parameterList$project_id))){
dir.create(paste0(Obio@parameterList$localWorkDir,Obio@parameterList$project_id))
}
if (!dir.exists(Obio@parameterList$reportFigDir)){
dir.create(Obio@parameterList$reportFigDir)
}
figureCount <- 1
## Load R module load R/3.5.1-foss-2018b ##
#setwd(Obio@parameterList$localWorkDir)
if (Obio@parameterList$host == "10.27.241.234"){
urlString <- "biologic.thecrick.org"
} else {
urlString <- "biologic.crick.ac.uk"
}###############################################################################
## Compile reference gene list ##
## Set gene reference list
Obio@dataTableList[["referenceList"]] <- list()
if (is.null(Obio@parameterList$HmDisplayCatsFromDb)){
Obio@parameterList$HmDisplayCatsFromDb <- list("TFs" = "ag_lab_categories__10")
}
## Add transcription factors that happen to be cluster markers ##
catList <- Obio@parameterList$HmDisplayCatsFromDb
if (Obio@parameterList$geneIDcolumn != "mgi_symbol" & Obio@parameterList$geneIDcolumn != "hgnc_symbol") {
queryGS <- "hgnc_symbol"
} else {
queryGS <- Obio@parameterList$geneIDcolumn
}
for (i in 1:length(catList)){
tempVec <- retrieve.gene.category.from.db(
cat_id = catList[[i]],
password = db.pwd,
gene.symbol = queryGS,
user = Obio@parameterList$db.user,
host = Obio@parameterList$host
)
Obio@dataTableList$referenceList[[names(catList)[i]]] <- tempVec
}
## Transcription Factors in variable genes ##
#tempVec <- tempVec[tempVec %in% Obio@dataTableList$dfGeneralMarkers$gene]Create Sample List
###############################################################################
## Create sample list filtered on MT and norm_counts_RNA ##
SampleList <- createNormSampleList(
obj = Obio,
reduce = NULL
)
## Done ##
###############################################################################
###############################################################################
## Add doublet annotation, if present, to meta data ##
pos <- grep("DF_resultlist", names(Obio@dataTableList))
if (length(pos) > 0){
sampleNames <- names(SampleList)
for (i in 1:length(SampleList)){
dfAdd <- Obio@dataTableList[["DF_resultlist"]][[sampleNames[i]]]
row.names(dfAdd) <- gsub("-1", "",row.names(dfAdd))
dfAdd <- dfAdd[row.names(dfAdd) %in% row.names(SampleList[[i]]@meta.data),]
SampleList[[i]] <- addDf2seuratMetaData(
obj = SampleList[[i]],
dfAdd = dfAdd
)
}
}
## Done ##
###############################################################################
###############################################################################
## Integrate Datasets ##
if (length(SampleList) > 1){
if (Obio@parameterList$scIntegrationMethod == "SCT"){
if (length(grep("scNintegrationFeatures", names(Obio@parameterList))) == 0){
Obio@parameterList$scNintegrationFeatures = 3000
}
library(future)
options(future.globals.maxSize = 7000 * 1024^2)
plan("multiprocess", workers = 30)
sample.features <- SelectIntegrationFeatures(
object.list = SampleList,
nfeatures = Obio@parameterList$scNintegrationFeatures
)
SampleList <- PrepSCTIntegration(
object.list = SampleList,
anchor.features = sample.features,
verbose = FALSE
)
sampleAnchors <- FindIntegrationAnchors(
object.list = SampleList,
normalization.method = "SCT",
anchor.features = sample.features,
verbose = FALSE
)
OsC <- IntegrateData(
anchorset = sampleAnchors,
normalization.method = "SCT",
verbose = FALSE
)
detach("package:future", unload=TRUE)
} else {
sampleAnchors <- FindIntegrationAnchors(
object.list = SampleList,
dims = 1:30
)
OsC <- IntegrateData(
#features.to.integrate = geneIntersectVec,
anchorset = sampleAnchors,
dims = 1:30
)
}
Obio@dataTableList$referenceList[["sampleAnchors"]] <- as.vector(sort(sampleAnchors@anchor.features))
} else {
OsC <- SampleList[[1]]
}
Idents(OsC) <- factor(Idents(OsC), levels = names(Obio@sampleDetailList))
dotsize = 1
if (nrow(OsC@meta.data) > 10000){
dotsize = 0.75
} else if (nrow(OsC@meta.data) > 20000){
dotsize = 0.5
} else if (nrow(OsC@meta.data) > 50000){
dotsize = 0.25
}
## Done integration ##
###############################################################################
###############################################################################
## Add additional annotation data ##
# ###############################################################################
# ## Load extra annotation data ##
# FN <- "/camp/stp/babs/working/boeings/Projects/swantonc/dhruva.biswas/335_CSL_DB_scRNA_seq_lung_tumor_microenvironment_EMTAB6653/basedata/MetaData.txt"
#
# dfExtra <- read.delim(
# FN,
# header = T,
# sep = "\t",
# stringsAsFactors = F
# )
#
# dfExtra <- dfExtra[dfExtra$cell %in% row.names(OsC@meta.data),]
#
# names(dfExtra) <- gsub("^cluster$", "Article_Cluster", names(dfExtra))
# names(dfExtra) <- gsub("^CellType$", "Article_Cell_Type", names(dfExtra))
# names(dfExtra) <- gsub("^CellType$", "Article_Cell_Type", names(dfExtra))
# names(dfExtra) <- gsub("^CellType$", "Article_Cell_Type", names(dfExtra))
#
# dfExtra[["Patient"]] <- sapply(dfExtra$Patient_piece, function(x) unlist(strsplit(x, "_"))[1])
# dfExtra$Patient <- paste0("P", dfExtra$Patient)
# dfExtra[["Region"]] <- sapply(dfExtra$Patient_piece, function(x) unlist(strsplit(x, "_"))[2])
# dfExtra$CellFromTumor <- as.character(dfExtra$CellFromTumor)
# dfExtra$CellFromTumor[dfExtra$CellFromTumor == "TRUE"] <- "Tumor"
# dfExtra$CellFromTumor[dfExtra$CellFromTumor == "FALSE"] <- "Non-Tumor"
#
#
#
# row.names(dfExtra) <- dfExtra$cell
#
# selVec <- c("CellFromTumor", "Patient", "Region")
# dfAdd <- dfExtra[,selVec]
#
# OsC <- addDf2seuratMetaData(
# obj = OsC,
# dfAdd = dfAdd
# )
#
# ## Done extra annotation ##
# ###############################################################################Result Figures
Quality Control Plots: Variable Features in the combined and all individual dataset
plotList <- list()
chnkVec <- as.vector(NULL, mode = "character")
## First make variation plot for integrated samples, than for all individual samples separately
tag <- "Integrated_Samples"
DefaultAssay(OsC) <- "RNA"
OsC <- FindVariableFeatures(
object = OsC,
selection.method = 'vst',
nfeatures = 2000
)
# Identify the 10 most highly variable genes
label2000 <- paste0("integrated", "_", "top2000var")
Obio@dataTableList$referenceList[[ label2000]]<- head(
x = VariableFeatures(object = OsC),
2000
)
label30 <- paste0("integrated", "_", "top30var")
Obio@dataTableList$referenceList[[ label30]]<- head(
x = VariableFeatures(object = OsC),
30
)
## slot for variable features OsC@assays$RNA@var.features
dfVar <- OsC@assays$RNA@meta.features
names(dfVar) <- gsub("vst.", "",names(dfVar))
dfVar[["gene"]] <- row.names(dfVar)
OsC@meta.data[["all"]] <- "all"
Idents(OsC) <- "all"
cluster.averages <- AverageExpression(
OsC,
return.seurat = TRUE
)
Idents(OsC) <- "sampleID"
dfAvgExpr <- data.frame(cluster.averages[["RNA"]]@data)
dfAvgExpr[["gene"]] <- row.names(dfAvgExpr)
names(dfAvgExpr)[1] <- "Avg.Expression"
dfVar <- merge(
dfVar,
dfAvgExpr,
by.x = "gene",
by.y = "gene"
)
dfVar[["Type"]] <- "Standard"
dfVar[dfVar$gene %in% OsC@assays$RNA@var.features, "Type"] <- "Most Variable"
dfVar[["text"]] <- ""
dfVar[dfVar$gene %in% as.vector(Obio@dataTableList$referenceList[[label30]]), "text"] <- Obio@dataTableList$referenceList[label30]
#dotsize <- 0.5
library(ggrepel)
plotList[[tag]] <- ggplot(
data = dfVar,
aes(
x=Avg.Expression,
y=variance.standardized, label = text, color = Type
)
) + geom_point( shape=16, size = dotsize
) + xlab("Average Expression") + ylab("Variance Standarized") + theme(
axis.text.y = element_text(size=8),
axis.text.x = element_text(size=8),
axis.title.y = element_text(size=8),
axis.title.x = element_text(size=8),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title = element_text(hjust = 0.5, size = 12)
) + ggtitle(paste0("Variance vs. Expression in the Overall Experiment")
) + scale_color_manual(values=c("#FF0000", "#000000")
) + geom_text_repel()
###########################################################################
## Save plot to file ##
FNbase <- paste0("variation.integrated.samples.", VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotList[[tag]])
dev.off()
## ##
###########################################################################
link <- paste0('< a href="https://',urlString,'/',Obio@parameterList$project_id,'/scatterplot?x_axis=',paste0(tag, '_AvgExpr'),'&y_axis=',paste0(tag, '_var_std" target = "_blank">here</a>'))
figCap <- paste0(
'**Figure ',
figureCount,
':** Variance versus averaged gene expression for overall sample.',
names(SampleList)[i],
'. ',
'Download a pdf of this figure <a href="',FNrel,'" target="_blank">here</a>. ',
'An interactive version of this figure can be found ', link, '. '
)
NewChnk <- paste0(
"#### ",tag,
"\n```{r varplot_",tag,", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",figCap,"'}\n",
"\n",
"\n print(plotList[['",tag,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
## Histogram Part C done ##
###########################################################################
chnkVec <- c(
chnkVec,
NewChnk
)
figureCount <- figureCount + 1
print(" All variation done.")
## Now the individual samples ##
xmax <- ceiling(max(dfVar$Avg.Expression))
ymax <- ceiling(max(dfVar$variance.standardized))
dfVarRes <- unique(
dfVar[,c("gene", "Avg.Expression", "variance.standardized")]
)
names(dfVarRes) <- gsub(
"Avg.Expression", paste0(tag, "_AvgExpr"), names(dfVarRes)
)
names(dfVarRes) <- gsub(
"variance.standardized", paste0(tag, "_var_std"), names(dfVarRes)
)
for (i in 1:length(SampleList)){
tag <- paste0("Ind_var_", names(SampleList)[i])
DefaultAssay(SampleList[[i]]) <- "RNA"
SampleList[[i]] <- FindVariableFeatures(object = SampleList[[i]])
SampleList[[i]] <- FindVariableFeatures(
object = SampleList[[i]],
selection.method = 'vst',
nfeatures = 2000
)
# Identify the 10 most highly variable genes
label2000 <- paste0(names(SampleList)[i], "_", "top2000var")
Obio@dataTableList$referenceList[[ label2000]]<- head(
x = VariableFeatures(object = SampleList[[i]]),
2000
)
label30 <- paste0(names(SampleList)[i], "_", "top30var")
Obio@dataTableList$referenceList[[ label30]]<- head(
x = VariableFeatures(object = SampleList[[i]]),
30
)
## slot for variable features OsC@assays$RNA@var.features
dfVar <- SampleList[[i]]@assays$RNA@meta.features
names(dfVar) <- gsub("vst.", "",names(dfVar))
dfVar[["gene"]] <- row.names(dfVar)
cluster.averages <- AverageExpression(
SampleList[[i]],
return.seurat = TRUE
)
dfAvgExpr <- data.frame(cluster.averages[["RNA"]]@data)
dfAvgExpr[["gene"]] <- row.names(dfAvgExpr)
names(dfAvgExpr)[1] <- "Avg.Expression"
dfVar <- merge(
dfVar,
dfAvgExpr,
by.x = "gene",
by.y = "gene"
)
dfVar[["Type"]] <- "Standard"
dfVar[dfVar$gene %in% SampleList[[i]]@assays$RNA@var.features, "Type"] <- "Most Variable"
dfVar[["text"]] <- ""
dfVar[dfVar$gene %in% as.vector(Obio@dataTableList$referenceList[[label30]]), "text"] <- Obio@dataTableList$referenceList[label30]
# dotsize <- 0.5
library(ggrepel)
plotList[[tag]] <- ggplot(
data = dfVar,
aes(
x=Avg.Expression,
y=variance.standardized, label = text, color = Type)
) + geom_point( shape=16, size = dotsize
) + xlab("Average Expression") + ylab("Variance Standarized") + theme(
axis.text.y = element_text(size=8),
axis.text.x = element_text(size=8),
axis.title.y = element_text(size=8),
axis.title.x = element_text(size=8),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title = element_text(hjust = 0.5, size = 12)
)+ ggtitle(paste0("Individual variance in the ", names(SampleList)[i], " sample.")
) + scale_color_manual(values=c("#FF0000", "#000000")
) + geom_text_repel(
) + xlim(0, xmax) + ylim(0, ymax)
###########################################################################
## Save plot to file ##
FNbase <- paste0("Individual.var.features",names(SampleList)[i], VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotList[[tag]])
dev.off()
## ##
###########################################################################
link <- paste0('<a href="https://',urlString,'/',Obio@parameterList$project_id,'/scatterplot?x_axis=',paste0(tag, '_AvgExpr'),'&y_axis=',paste0(tag, '_var_std" target="_blank">here</a>'))
figCap <- paste0(
'**Figure ',
figureCount,
':** Variance versus averaged gene expression for sample ',
names(SampleList)[i],
'. ',
'Download a pdf of this figure <a href="',FNrel,'" target="_blank">here</a>. ',
'An interactive version of this figure can be found ', link, '. '
)
NewChnk <- paste0(
"#### ",tag,
"\n```{r varplot_",tag,", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",figCap,"'}\n",
"\n",
"\n print(plotList[['",tag,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
## Histogram Part C done ##
###########################################################################
chnkVec <- c(
chnkVec,
NewChnk
)
figureCount <- figureCount + 1
print(paste0(names(SampleList)[i], " individual variation done."))
## Add to result table ##
dfVarTemp <- unique(
dfVar[,c("gene", "Avg.Expression", "variance.standardized")]
)
names(dfVarTemp) <- gsub(
"Avg.Expression", paste0(tag, "_AvgExpr"), names(dfVarTemp)
)
names(dfVarTemp) <- gsub(
"variance.standardized", paste0(tag, "_var_std"), names(dfVarTemp)
)
dfVarRes <- merge(
dfVarRes,
dfVarTemp,
by.x = "gene",
by.y = "gene",
all = TRUE
)
dfVarRes[is.na(dfVarRes)] <- 0
}
## Make sure summary goes first ##
Obio@dataTableList[["dfVariation"]] <- dfVarRes## plot list will be integrated in full figure ##
cat(paste(knit(text = chnkVec, quiet = T), collapse = '\n'))Integrated_Samples
Figure 1: Variance versus averaged gene expression for overall sample.S. Download a pdf of this figure here. An interactive version of this figure can be found < a href=“https://biologic.crick.ac.uk/vpl362A/scatterplot?x_axis=Integrated_Samples_AvgExpr&y_axis=Integrated_Samples_var_std” target = "_blank">here.
Ind_var_S
Figure 2: Variance versus averaged gene expression for sample S. Download a pdf of this figure here. An interactive version of this figure can be found here.
###############################################################################
## Perform integrated analysis ##
if (length(Obio@sampleDetailList) > 1){
DefaultAssay(OsC) <- "integrated"
} else {
Obio@parameterList$singleCellClusterString <- gsub("integrated", "RNA", Obio@parameterList$singleCellClusterString)
}
# Run the standard workflow for visualization and clustering
## This will scale on the most variable features only
OsC <- ScaleData(OsC, verbose = FALSE)
OsC <- RunPCA(
OsC,
npcs = Obio@parameterList$singleCellSeuratNpcs4PCA, verbose = FALSE
)
# t-SNE and Clustering
## Add PCA clusters to data collection ##
OsC <- RunUMAP(OsC, reduction = "pca", dims = 1:20)
OsC <- RunTSNE(OsC, reduction = "pca", dims = 1:20)
OsC <- FindNeighbors(OsC, reduction = "pca", dims = 1:20)
OsC <- FindClusters(OsC, resolution = Obio@parameterList$singleCellClusterParameter)
## Rational: Run PCA on variable features, then scale data for heatmaps and other applications
if (Obio@parameterList$scIntegrationMethod != "RNA"){
DefaultAssay(OsC) <- Obio@parameterList$scIntegrationMethod
allGenes <- rownames(x = OsC@assays[[Obio@parameterList$scIntegrationMethod]])
OsC <- ScaleData(OsC, verbose = FALSE, features=allGenes)
}
DefaultAssay(OsC) <- "RNA"
allGenes <- rownames(x = OsC@assays$RNA)
OsC <- ScaleData(OsC, verbose = FALSE, features=allGenes)Dimensionality Reduction Plots by Sample
If you wish to get a bit of background on tSNE dimensionality reduction, take a look at this youtube video by Josh Starmer from the University of North Carolina.
If you wish to get a bit of background on UMAP (and other) dimensionality reduction algorithms, take a look at this youtube video by recaping a lecture at the PyData 2018 conference.
reductionVec <- c("umap", "tsne")
plotList <- list()
chnkVec <- as.vector(NULL, mode = "character")
## Merge UMAP Coordinates and cluster assignment ##
for (i in 1:length(reductionVec)){
tag <- paste0("Dimplot_by_sample_", i)
plotList[[tag]] <- DimPlot(OsC, reduction = reductionVec[i], group.by = "sampleID")
## Save to file ##
FNbase <- paste0("dimplot.by.sample.", reductionVec[i],".", VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotList[[tag]])
dev.off()
if (reductionVec[i] == "tsne"){
link <- paste0('<a href="https://',urlString,'/',Obio@parameterList$project_id,'/pca?x_axis=tSNE_1&y_axis=tSNE_2" target="_blank">here</a>')
} else {
link <- paste0('<a href="https://',urlString,'/',Obio@parameterList$project_id,'/pca?x_axis=UMAP_1&y_axis=UMAP_2" target="_blank">here</a>')
}
figLegend <- paste0(
'**Figure ',
figureCount,
':** ',
reductionVec[i],
' plot depicting all samples. Download a pdf of this figure <a href="',FNrel,'" target="_blank">here</a>.',
'An interactive version of this figure can be found ', link, '. '
)
NewChnk <- paste0(
"#### ",reductionVec[i],
"\n```{r Dimplot_by_sample_",
i,", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",
figLegend,"'}\n",
"\n",
"\n print(plotList[['",tag,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
chnkVec <- c(
chnkVec,
NewChnk
)
}
## Done integraed analysis ##
###############################################################################cat(paste(knit(text = chnkVec, quiet = T), collapse = '\n'))umap
Calculate Average Expression Per Cluster
## Ensure that hmIdent exists ##
## Want to average by cluster and subset by sampleID ##
###############################################################################
## Calculataing average expression across all samples ##
OsC@meta.data[["all"]] <- "all"
Idents(OsC) <- "all"
cluster.averages <- AverageExpression(OsC, return.seurat = TRUE)
## Retrieved Scaled data ##
dfAvgExpr <- data.frame(cluster.averages[["RNA"]]@data)
#dfAvgExpr <- dfAvgExpr[,sort(names(dfAvgExpr))]
dfAvgExpr[["gene"]] <- row.names(dfAvgExpr)
dfAvgScaledData <- data.frame(cluster.averages[["RNA"]]@scale.data)
#dfAvgScaledData <- dfAvgScaledData[,sort(names(dfAvgScaledData))]
dfAvgScaledData[["gene"]] <- row.names(dfAvgScaledData)
Obio@dataTableList[["dfAvglg10ExprAll"]] <- dfAvgExpr
Obio@dataTableList[["dfAvglg10ExprAllScaled"]] <- dfAvgScaledData
## Done ##
###############################################################################
OsC@meta.data[["clustIdent"]] <- paste0(
"C_", OsC@meta.data[,Obio@parameterList$singleCellClusterString]
)
Idents(OsC) <- "clustIdent"
if (length(grep("sampleID", names(OsC@meta.data))) > 0){
cluster.averages <- AverageExpression(OsC, return.seurat = TRUE, add.ident = "sampleID")
} else {
cluster.averages <- AverageExpression(OsC, return.seurat = TRUE)
}
Idents(OsC) <- "sampleID"
## Retrieved Scaled data ##
dfAvgExpr <- data.frame(cluster.averages[["RNA"]]@data)
dfAvgExpr <- dfAvgExpr[,sort(names(dfAvgExpr))]
dfAvgExpr[["gene"]] <- row.names(dfAvgExpr)
dfAvgScaledData <- data.frame(cluster.averages[["RNA"]]@scale.data)
dfAvgScaledData <- dfAvgScaledData[,sort(names(dfAvgScaledData))]
dfAvgScaledData[["gene"]] <- row.names(dfAvgScaledData)
Obio@dataTableList[["dfAvglg10ExprByClusterBySample"]] <- dfAvgExpr
Obio@dataTableList[["dfAvglg10ExprByClusterBySampleScaled"]] <- dfAvgScaledData
###############################################################################
###############################################################################
## Average by Cluster ##
Idents(OsC) <- "clustIdent"
Idents(OsC) <- factor(Idents(OsC), levels = sort(levels(OsC)))
cluster.averages <- AverageExpression(OsC, return.seurat = TRUE)
dfAvgExpr <- data.frame(cluster.averages[["RNA"]]@data)
dfAvgExpr <- dfAvgExpr[,sort(names(dfAvgExpr))]
dfAvgExpr[["gene"]] <- row.names(dfAvgExpr)
dfAvgScaledData <- data.frame(cluster.averages[["RNA"]]@scale.data)
dfAvgScaledData <- dfAvgScaledData[,sort(names(dfAvgScaledData))]
dfAvgScaledData[["gene"]] <- row.names(dfAvgScaledData)
Obio@dataTableList[["dfAvglg10ExprPerCluster"]] <- dfAvgExpr
Obio@dataTableList[["dfAvglg10ExprPerClusterScaled"]] <- dfAvgScaledData
## Done Average by Cluster ##
###############################################################################
###############################################################################
## Average gene expression by sample ##
Idents(OsC) <- "sampleID"
cluster.averages <- AverageExpression(OsC, return.seurat = TRUE)
## Retrieved Scaled data ##
dfAvgExpr <- data.frame(cluster.averages[["RNA"]]@data)
selVec <- names(dfAvgExpr)
dfAvgExpr[["gene"]] <- row.names(dfAvgExpr)
selVec <- c("gene", selVec)
dfAvgExpr <- dfAvgExpr[,selVec]
dfAvgScaledData <- data.frame(cluster.averages[["RNA"]]@scale.data)
selVec <- names(dfAvgScaledData)
dfAvgScaledData[["gene"]] <- row.names(dfAvgScaledData)
selVec <- c("gene", selVec)
dfAvgScaledData <- dfAvgScaledData[,selVec]
Obio@dataTableList[["dfAvglg10ExprBySample"]] <- dfAvgExpr
Obio@dataTableList[["dfAvglg10ExprBySampleScaled"]] <- dfAvgScaledData
## Done average gene expression by sample ##
###############################################################################Dimensionality Reduction Plots by Cluster
reductionVec <- c("umap", "tsne")
plotList <- list()
chnkVec <- as.vector(NULL, mode = "character")
###############################################################################
## First UMAP all samples together ##
tag <- paste0("UMAP_All_Samples")
dfPlot <- OsC@meta.data
pos <- grep("included", names(dfPlot))
if (length(pos) == 0){
dfPlot[["included"]] <- "+"
}
dfPlot[["cellID"]] <- row.names(dfPlot)
dfPlot$UMAP_1 <- NULL
dfPlot$UMAP_2 <- NULL
## Get UMAP coordinates ##
coord <- data.frame(OsC@reductions$umap@cell.embeddings)
coord[["cellID"]] <- row.names(coord)
coord <-coord[coord$cellID %in% dfPlot$cellID, ]
dfPlot <- merge(dfPlot, coord, by.x = "cellID", by.y="cellID", all=T)
dfPlot[is.na(dfPlot)] <- 0
dfPlot <- dfPlot[dfPlot$UMAP_1 != 0 & dfPlot$UMAP_2 != 0,]
## Add cluster colors ##
dfPlot[["Cluster"]] <- paste0("C", dfPlot$seurat_clusters)
clusterVec <- as.vector(paste0("C", unique(sort(dfPlot$seurat_clusters))))
maxX <- 1.1*max(dfPlot$UMAP_1, na.rm = T)
minX <- 1.1*min(dfPlot$UMAP_1, na.rm = T)
maxY <- 1.1*max(dfPlot$UMAP_2, na.rm = T)
minY <- 1.1*min(dfPlot$UMAP_2, na.rm = T)
library(scales)
clusterCols = hue_pal()(length(clusterVec))
dfPlot$Cluster <- factor(dfPlot$Cluster, levels = clusterVec)
# dotsize = 1
# if (nrow(dfPlot) > 10000){
# dotsize = 0.75
# } else if (nrow(dfPlot) > 20000){
# dotsize = 0.5
# } else if (nrow(dfPlot) > 50000){
# dotsize = 0.25
# }
plotList[[tag]] <- ggplot(data=dfPlot[dfPlot$included == "+",], aes(UMAP_1, UMAP_2, color=Cluster)
) + geom_point( shape=16, size = as.numeric(dotsize)
) + xlab("UMAP1") + ylab("UMAP2") + theme(
axis.text.y = element_text(size=8),
axis.text.x = element_text(size=8),
axis.title.y = element_text(size=8),
axis.title.x = element_text(size=8),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title = element_text(hjust = 0.5, size = 12),
legend.title = element_blank()
) + ggtitle(paste0("Sample: ", tag)
) + xlim(minX, maxX) + ylim(minY, maxY
) + coord_fixed(ratio=1
) + theme_bw()
if (length(unique(dfPlot$Cluster)) > 15){
plotList[[tag]] <- plotList[[tag]] + theme(legend.position = "none")
}
FNbase <- paste0(tag, VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotList[[tag]])
dev.off()
link <- paste0('<a href="https://',urlString,'/',Obio@parameterList$project_id,'/pca?x_axis=UMAP_1&y_axis=UMAP_2" target="_blank">here</a>')
figLegend <- paste0(
'**Figure ',
figureCount,
':** ',
' UMAP showing all cells from all samples together. Download a pdf of this figure <a href="',FNrel,'" target="_blank">here</a>.',
'An interactive version of this figure can be found ', link, '. '
)
figureCount <- figureCount + 1
NewChnk <- paste0(
"#### ", tag,
"\n```{r SL_UMAP_",
tag,", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",
figLegend,"'}\n",
"\n",
"\n print(plotList[['",tag,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
chnkVec <- c(
chnkVec,
NewChnk
)
## Done first umap all samples ##
###############################################################################
###############################################################################
## First tsne all samples together ##
tag <- paste0("tSNE_All_Samples")
dfPlot <- OsC@meta.data
pos <- grep("included", names(dfPlot))
if (length(pos) == 0){
dfPlot[["included"]] <- "+"
}
dfPlot[["cellID"]] <- row.names(dfPlot)
dfPlot$tSNE_1 <- NULL
dfPlot$tSNE_2 <- NULL
## Get tsNE coordinates ##
coord <- data.frame(OsC@reductions$tsne@cell.embeddings)
coord[["cellID"]] <- row.names(coord)
coord <-coord[coord$cellID %in% dfPlot$cellID, ]
dfPlot <- merge(dfPlot, coord, by.x = "cellID", by.y="cellID", all=T)
dfPlot[is.na(dfPlot)] <- 0
dfPlot <- dfPlot[dfPlot$tSNE_1 != 0 & dfPlot$tSNE_2 != 0,]
## Add cluster colors ##
dfPlot[["Cluster"]] <- paste0("C", dfPlot$seurat_clusters)
clusterVec <- as.vector(paste0("C", unique(sort(dfPlot$seurat_clusters))))
maxX <- 1.1*max(dfPlot$tSNE_1, na.rm = T)
minX <- 1.1*min(dfPlot$tSNE_1, na.rm = T)
maxY <- 1.1*max(dfPlot$tSNE_2, na.rm = T)
minY <- 1.1*min(dfPlot$tSNE_2, na.rm = T)
library(scales)
clusterCols = hue_pal()(length(clusterVec))
dfPlot$Cluster <- factor(dfPlot$Cluster, levels = clusterVec)
# dotsize = 1.5
# if (nrow(dfPlot) > 10000){
# dotsize = 0.75
# } else if (nrow(dfPlot) > 50000){
# dotsize = 0.5
# } else {
# dotsize = 0.25
# }
plotList[[tag]] <- ggplot(data=dfPlot[dfPlot$included == "+",], aes(tSNE_1, tSNE_2, color=Cluster)
) + geom_point( shape=16, size = as.numeric(dotsize)
) + xlab("tSNE1") + ylab("tSNE2") + theme(
axis.text.y = element_text(size=8),
axis.text.x = element_text(size=8),
axis.title.y = element_text(size=8),
axis.title.x = element_text(size=8),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title = element_text(hjust = 0.5, size = 12),
legend.title = element_blank()
) + ggtitle(paste0("Sample: ", tag)
) + xlim(minX, maxX) + ylim(minY, maxY
) + coord_fixed(ratio=1
) + theme_bw()
if (length(unique(dfPlot$Cluster)) > 15){
plotList[[tag]] <- plotList[[tag]] + theme(legend.position = "none")
}
FNbase <- paste0(tag, VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotList[[tag]])
dev.off()
link <- paste0('<a href="https://',urlString,'/',Obio@parameterList$project_id,'/pca?x_axis=tSNE_1&y_axis=tSNE_2" target="_blank">here</a>')
figLegend <- paste0(
'**Figure ',
figureCount,
':** ',
' tSNE showing all cells from all samples together. Download a pdf of this figure <a href="',FNrel,'" target="_blank">here</a>.',
'An interactive version of this figure can be found ', link, '. '
)
figureCount <- figureCount + 1
NewChnk <- paste0(
"#### ", tag,
"\n```{r SL_tSNE_",
tag,", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",
figLegend,"'}\n",
"\n",
"\n print(plotList[['",tag,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
chnkVec <- c(
chnkVec,
NewChnk
)
## Done first tsne all samples ##
###############################################################################
###############################################################################
## Make one UMAP plot per sample ##
sampleVec <- sort(unique(OsC@meta.data$sampleID))
dfPlot <- OsC@meta.data
pos <- grep("included", names(dfPlot))
if (length(pos) == 0){
dfPlot[["included"]] <- "+"
}
dfPlot[["cellID"]] <- row.names(dfPlot)
## Get UMAP coordinates ##
coord <- data.frame(OsC@reductions$umap@cell.embeddings)
coord[["cellID"]] <- row.names(coord)
coord <-coord[coord$cellID %in% dfPlot$cellID, ]
dfPlot$UMAP_1 <- NULL
dfPlot$UMAP_2 <- NULL
dfPlot <- merge(dfPlot, coord, by.x = "cellID", by.y="cellID", all=T)
dfPlot[is.na(dfPlot)] <- 0
dfPlot <- dfPlot[dfPlot$UMAP_1 != 0 & dfPlot$UMAP_2 != 0,]
## Add cluster colors ##
dfPlot[["Cluster"]] <- paste0("C", dfPlot$seurat_clusters)
clusterVec <- as.vector(paste0("C", unique(sort(dfPlot$seurat_clusters))))
library(scales)
clusterCols = hue_pal()(length(clusterVec))
dfPlot$Cluster <- factor(dfPlot$Cluster, levels = clusterVec)
maxX <- 1.1*max(dfPlot$UMAP_1, na.rm = T)
minX <- 1.1*min(dfPlot$UMAP_1, na.rm = T)
maxY <- 1.1*max(dfPlot$UMAP_2, na.rm = T)
minY <- 1.1*min(dfPlot$UMAP_2, na.rm = T)
for (i in 1:length(sampleVec)){
tag <- paste0("UMAP_plot_by_", sampleVec[i])
dfPlotSel <- dfPlot[dfPlot$sampleID == sampleVec[i], ]
plotList[[tag]] <- ggplot(data=dfPlotSel[dfPlot$included == "+",], aes(UMAP_1, UMAP_2, color=Cluster)
) + geom_point( shape=16, size = as.numeric(dotsize)
) + xlab("UMAP1") + ylab("UMAP2") + theme(
axis.text.y = element_text(size=8),
axis.text.x = element_text(size=8),
axis.title.y = element_text(size=8),
axis.title.x = element_text(size=8),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title = element_text(hjust = 0.5, size = 12),
legend.title = element_blank()
) + ggtitle(paste0("Sample: ", tag)
) + xlim(minX, maxX) + ylim(minY, maxY
) + coord_fixed(ratio=1
) + theme_bw()
if (length(unique(dfPlot$Cluster)) > 15){
plotList[[tag]] <- plotList[[tag]] + theme(legend.position = "none")
}
FNbase <- paste0(tag, VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotList[[tag]])
dev.off()
figLegend <- paste0(
'**Figure ',
figureCount,
':** ',
' Sample-level UMAPs. Download a pdf of this figure <a href="',FNrel,'" target="_blank">here</a>.'
)
figureCount <- figureCount + 1
NewChnk <- paste0(
paste("#### ", tag),
"\n```{r SL_UMAP_",
tag,", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",
figLegend,"'}\n",
"\n",
"\n print(plotList[['",tag,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
chnkVec <- c(
chnkVec,
NewChnk
)
}
## Done making one umap plot per sample ##
###############################################################################
###############################################################################
## Add cluster dendrogram by sample ##
if (length(unique(OsC@meta.data$sampleID)) > 3){
library(ggtree)
Idents(OsC) <- "sampleID"
OsC <- BuildClusterTree(OsC)
tag <- paste0("Sample_Dendrogram")
OsC@tools$BuildClusterTree$tip.label <- paste0("C", OsC@tools$BuildClusterTree$tip.label)
plotList[[tag]] <- ggplot(OsC@tools$BuildClusterTree
) + geom_tree(
) + theme_tree(
) + geom_tiplab(
) + labs(title=tag
) + theme(
panel.border = element_rect(colour = "black", fill=NA, size=1),
axis.title.x=element_blank(),
plot.title = element_text(hjust = 0.5, size = 12)
) + xlim(0,1.2*max(OsC@tools$BuildClusterTree[[2]]))
## Save to file ##
FNbase <- paste0(tag,".", VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotList[[tag]])
dev.off()
figLegend <- paste0(
'**Figure ',
figureCount,
':** ',
' Clusterplot dendrogram by sample ID. ','A pdf of this figure can be downloaded <a href="',FNrel,'" target="_blank">here</a>.'
)
NewChnk <- paste0(
"#### SampleID Dendrogram",
"\n```{r ", tag, "results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",
figLegend,"'}\n",
"\n",
"\n print(plotList[['",tag,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
chnkVec <- c(
chnkVec,
NewChnk
)
figureCount <- figureCount + 1
}
## Done by sample ##
###############################################################################
###############################################################################
## Add cluster dendrogram by cluster ##
library(ggtree)
Idents(OsC) <- "seurat_clusters"
OsC <- BuildClusterTree(OsC)
tag <- paste0("Cluster_Dendrogram")
OsC@tools$BuildClusterTree$tip.label <- paste0("C", OsC@tools$BuildClusterTree$tip.label)
plotList[[tag]] <- ggplot(OsC@tools$BuildClusterTree
) + geom_tree(
) + theme_tree(
) + geom_tiplab(
) + labs(title=tag
) + theme(
panel.border = element_rect(colour = "black", fill=NA, size=1),
axis.title.x=element_blank(),
plot.title = element_text(hjust = 0.5, size = 12)
) + xlim(0, 2*max(OsC@tools$BuildClusterTree[[2]]))
#+ xlim(-1,1.2*max(OsC@tools$BuildClusterTree$edge))
## Save to file ##
FNbase <- paste0(tag,".", VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotList[[tag]])
dev.off()
figLegend <- paste0(
'**Figure ',
figureCount,
':** ',
' Clusterplot dendrogram. ','A pdf of this figure can be downloaded <a href="',FNrel,'" target="_blank">here</a>.'
)
NewChnk <- paste0(
"#### Cluster Dendrogram",
"\n```{r ", tag, "results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",
figLegend,"'}\n",
"\n",
"\n print(plotList[['",tag,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
chnkVec <- c(
chnkVec,
NewChnk
)
figureCount <- figureCount + 1
## Done integraed analysis ##
###############################################################################
###############################################################################
## Find all markers ##
DefaultAssay(OsC) <- "RNA"
Idents(OsC) <- Obio@parameterList$singleCellClusterString
lgFCthreshold <- 0.25
dfGeneralMarkers <- FindAllMarkers(
object = OsC,
only.pos = FALSE,
min.pct = 0.1,
logfc.threshold = lgFCthreshold,
test.use = "roc",
assay = "RNA",
slot = "data"
)
if (nrow(dfGeneralMarkers) == 0){
lgFCthreshold <- 0.00
dfGeneralMarkers <- FindAllMarkers(
object = OsC,
only.pos = FALSE,
min.pct = 0.01,
logfc.threshold = lgFCthreshold,
test.use = "roc",
assay = "RNA",
slot = "data"
)
}
if (nrow(dfGeneralMarkers) > 0){
dfGeneralMarkers[["direction"]] <- ""
dfGeneralMarkers[dfGeneralMarkers$avg_diff >= 0, "direction"] <- "positive"
dfGeneralMarkers[dfGeneralMarkers$avg_diff < 0, "direction"] <- "negative"
Obio@dataTableList[["dfGeneralMarkers"]] <- dfGeneralMarkers
dfGeneralMarkersFilt <- dfGeneralMarkers[dfGeneralMarkers$avg_diff > lgFCthreshold | dfGeneralMarkers$avg_diff < -lgFCthreshold,]
Obio@dataTableList[["dfGeneralMarkersFilt"]] <- dfGeneralMarkersFilt
dfTop1 <- data.frame(dfGeneralMarkers %>% group_by(cluster) %>% top_n(1, avg_diff))
dfTop5 <- data.frame(dfGeneralMarkers %>% group_by(cluster) %>% top_n(5, avg_diff))
dfTop10 <- data.frame(dfGeneralMarkers %>% group_by(cluster) %>% top_n(10, avg_diff))
Obio@dataTableList[["dfGeneralMarkersTop10"]] <- dfTop10
Obio@dataTableList$referenceList[["Top10clusterMarkers"]] <- as.vector(
unique(
dfTop10$gene
)
)
} else {
Obio@dataTableList[["dfGeneralMarkers"]] <- NULL
}
#############################################################
###############################################################################
## Upload general markers ##
library(tidyverse)
dfDat <- dfGeneralMarkers
dfDat <- unique(dfGeneralMarkers[,c("gene", "cluster")])
if (Obio@parameterList$geneIDcolumn != "hgnc_symbol" & Obio@parameterList$geneIDcolumn != "mgi_symbol"){
refGeneIDcolumn <- "hgnc_symbol"
dfAnno <- Obio@dfGeneAnnotation
dfAnno <- unique(dfAnno[,c("hgnc_symbol",Obio@parameterList$geneIDcolumn )])
dfAnno <- dfAnno[dfAnno[,Obio@parameterList$geneIDcolumn] %in% dfDat[,"gene"],]
dfDat <- merge(
dfDat,
dfAnno,
by.x = "gene",
by.y = Obio@parameterList$geneIDcolumn
)
dfDat$gene <- NULL
names(dfDat) <- gsub("hgnc_symbol", "gene",names(dfDat))
} else {
refGeneIDcolumn <- Obio@parameterList$geneIDcolumn
}
dfDat <- dfDat %>%
group_by(cluster) %>%
mutate(rn = row_number()) %>%
ungroup %>%
pivot_wider(names_from = cluster, values_from = gene, values_fill = "") %>% select(-rn)
orderVec <- sort(names(dfDat))
dfDat <- dfDat[,orderVec]
names(dfDat) <- paste0(
Obio@parameterList$project_id,
"_Marker_Genes_Cluster_",
names(dfDat)
)
dfDat <- data.frame(dfDat)
## Insert description row ##
descriptionRow <- dfDat[1,]
descriptionRow[1,] <- t(names(dfDat))
descriptionRow[1,] <- paste0(descriptionRow[1,], " from Seurat FindAllMarkers.")
dfDat <- rbind(
descriptionRow,
dfDat
)
# dfCatRef <- names(dfCatRef)
# names(dfCatRef) <- names(dfDat)
#
# dfDat <- rbind(
# dfCatRef,
# dfDat
# )
#######################################################################
## Upload/update category by category ##
updatedCatIDs <- as.vector(NULL, mode = "character")
updatedCatNames <- as.vector(NULL, mode = "character")
for (i in 1:ncol(dfDat)){
cat.name <- names(dfDat)[i]
cat_type <- paste0("temp_cluster_marker_", Obio@parameterList$project_id)
cat.description.text <- as.vector(dfDat[1,i])
gene.vec <- as.vector(
dfDat[,i]
)[2:nrow(dfDat)]
gene.vec <- gene.vec[gene.vec != ""]
gene.vec <- sort(na.omit(gene.vec))
## Determine if cat exists ##
catID <- add.category.to.lab.reference.table.hs(
host = Obio@dbDetailList$host,
pwd = db.pwd,
user = Obio@dbDetailList$db.user,
cat.ref.db = Obio@dbDetailList$ref.cat.db,
cat.ref.db.table = Obio@parameterList$lab.categories.table,
gene.vector = gene.vec,
gene.id = refGeneIDcolumn, #options hgnc_symbol, mgi_symbol
mm.hs.conversion.file = paste0(hpc.mount, "Projects/reference_data/20160303.homologene.data.txt"),
cat_name = cat.name,
cat_type = cat_type,
data_source = paste0(Obio@parameterList$labname, " Lab"),
comments_1 = "",
comments_2 = "",
new.lab.category.table = FALSE,
cat.description.db = "internal_categories",
cat.description.db.table = "category_description",
cat.description.text = cat.description.text,
lab.name = Obio@parameterList$labname,
replaceExistingCatName = TRUE
)
updatedCatIDs <- c(
updatedCatIDs,
catID
)
updatedCatNames <- c(
updatedCatNames,
cat.name
)
} ## End dfDat loop
dfFAMplotIDs <- data.frame(
cat_id = updatedCatIDs,
cat_name = updatedCatNames
)
## Done uploading general markers ##
###############################################################################cat(paste(knit(text = chnkVec, quiet = T), collapse = '\n'))UMAP_All_Samples
tSNE_All_Samples
UMAP_plot_by_S
Figure 5: Sample-level UMAPs. Download a pdf of this figure here.
Cluster Dendrogram
Figure 6: Clusterplot dendrogram. A pdf of this figure can be downloaded here.
Check Accuracy UMAP Distances
The sleepwalk tool will provide the euclidean distances between individual cells. This will help you to determine clustering accuracy. Find more information here.library(sleepwalk)
sleepwalk(
OsC@reductions$umap@cell.embeddings,
OsC@reductions$pca@cell.embeddings,
saveToFile=paste(Obio@parameterList$outputDir,"sleepwalk.UMAP.html",sep='')
)
htmltools::includeHTML(paste(Obio@parameterList$outputDir,"sleepwalk.UMAP.html",sep=''))Highlight Potential Doublets
This analyis gives an indication of which cells MIGHT be doublets. Chris McGinnis DoubletFinder package is used to identify potential doublets.
## Add UMAP coordinates to Metadata ##
dfAdd <- data.frame(OsC@reductions$umap@cell.embeddings)
OsC <- addDf2seuratMetaData(
obj = OsC,
dfAdd = dfAdd
)
## Add tSNE coordinates to Metadata ##
dfAdd <- data.frame(OsC@reductions$tsne@cell.embeddings)
OsC <- addDf2seuratMetaData(
obj = OsC,
dfAdd = dfAdd
)
plotList <- list()
chnkVec <- as.vector(NULL, mode = "character")
dfTemp <- OsC@meta.data
pos <- grep("Doublet_plot", names(dfTemp))
if (length(pos) > 0){
## First make variation plot for integrated samples, than for all individual samples separately
tag <- "Doublet_plot"
dfTemp$DF_Classification <- factor(dfTemp$DF_Classification, levels = sort(unique(dfTemp$DF_Classification)))
#dotsize <- round(7500/nrow(dfTemp),1)
# dotsize <- 0.3
plotList[[tag]] <- ggplot(dfTemp, aes(UMAP_1, UMAP_2, color=DF_Classification)
)+ geom_point(
shape = 16,
size = as.numeric(dotsize)
) + xlab("UMAP1") + ylab("UMAP2") + theme(
axis.text.y = element_text(size=8),
axis.text.x = element_text(size=8),
axis.title.y = element_text(size=8),
axis.title.x = element_text(size=8),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title = element_text(hjust = 0.5, size = 12),
panel.background = element_rect(fill = "lightgrey")
) + ggtitle("Potential Doublets"
) + scale_color_manual(values=c("#000000","#FF0000")
) + theme_bw()
#+ xlim(minX, maxX) + ylim(minY, maxY)
## Save to file ##
FNbase <- paste0("DoubletFinderAll", VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotList[[tag]])
dev.off()
## Create R markdown chunk ##
figLegend <- paste0(
'**Figure ',
figureCount,
'**: Figure depicting the location of potential doublets in PCA components 1 and 2. Download a pdf of this figure <a href="',FNrel,'" target="_blank">here</a>. '
)
figureCount <- figureCount + 1
NewChnk <- paste0(
" #### Doublets All Timepoints",
"\n```{r ", tag, ", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",
figLegend,"'}\n",
"\n",
"\n print(plotList[['",tag,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
chnkVec <- c(
chnkVec,
NewChnk
)
}## plot list will be integrated in full figure ##
if (length(pos) > 0){
cat(paste(knit(text = chnkVec, quiet = T), collapse = '\n'))
}Cell Cycle Markers
###############################################################################
## Estimate cell cycle genes ##
exprFN <- paste0(hpc.mount, "Projects/reference_data/cell_cycle_vignette_files/nestorawa_forcellcycle_expressionMatrix.txt")
exp.mat <- read.table(file = exprFN, header = TRUE,
as.is = TRUE, row.names = 1)
# A list of cell cycle markers, from Tirosh et al, 2015, is loaded with Seurat. We can
# segregate this list into markers of G2/M phase and markers of S phase
s.genes <- cc.genes$s.genes
g2m.genes <- cc.genes$g2m.genes
print(paste0("Used as S-phase marker genes: ", sort(unique(paste(s.genes, collapse = ", ")))))
print(paste0("Used as G2M-phase marker genes: ", sort(unique(paste(g2m.genes, collapse = ", ")))))
# Create our Seurat object and complete the initalization steps
OsC <- CellCycleScoring(OsC, s.features = s.genes, g2m.features = g2m.genes, set.ident = TRUE)
reductionVec <- c("umap", "tsne")
plotList <- list()
chnkVec <- as.vector(NULL, mode = "character")
###############################################################################
## First UMAP all samples together ##
tag <- paste0("CellCyclePhase_All_Samples")
dfPlot <- OsC@meta.data
pos <- grep("included", names(dfPlot))
if (length(pos) == 0){
dfPlot[["included"]] <- "+"
}
dfPlot[["cellID"]] <- row.names(dfPlot)
dfPlot$UMAP_1 <- NULL
dfPlot$UMAP_2 <- NULL
## Get UMAP coordinates ##
coord <- data.frame(OsC@reductions$umap@cell.embeddings)
coord[["cellID"]] <- row.names(coord)
coord <-coord[coord$cellID %in% dfPlot$cellID, ]
dfPlot <- merge(dfPlot, coord, by.x = "cellID", by.y="cellID", all=T)
dfPlot[is.na(dfPlot)] <- 0
dfPlot <- dfPlot[dfPlot$UMAP_1 != 0 & dfPlot$UMAP_2 != 0,]
## Add cluster colors ##
#dfPlot[["Cluster"]] <- paste0("C", dfPlot$seurat_clusters)
#clusterVec <- as.vector(paste0("C", unique(sort(dfPlot$seurat_clusters))))
maxX <- 1.1*max(dfPlot$UMAP_1, na.rm = T)
minX <- 1.1*min(dfPlot$UMAP_1, na.rm = T)
maxY <- 1.1*max(dfPlot$UMAP_2, na.rm = T)
minY <- 1.1*min(dfPlot$UMAP_2, na.rm = T)
# library(scales)
# clusterCols = hue_pal()(length(clusterVec))
# dfPlot$Cluster <- factor(dfPlot$Cluster, levels = clusterVec)
# dotsize = 1
# if (nrow(dfPlot) > 10000){
# dotsize = 0.75
# } else if (nrow(dfPlot) > 20000){
# dotsize = 0.5
# } else if (nrow(dfPlot) > 50000){
# dotsize = 0.25
# }
plotList[[tag]] <- ggplot(data=dfPlot[dfPlot$included == "+",], aes(UMAP_1, UMAP_2, color=Phase)
) + geom_point( shape=16, size = as.numeric(dotsize)
) + xlab("UMAP1") + ylab("UMAP2") + theme(
axis.text.y = element_text(size=8),
axis.text.x = element_text(size=8),
axis.title.y = element_text(size=8),
axis.title.x = element_text(size=8),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title = element_text(hjust = 0.5, size = 12),
legend.title = element_blank()
) + ggtitle(paste0("Sample: ", tag)
) + xlim(minX, maxX) + ylim(minY, maxY
) + coord_fixed(ratio=1
) + theme_bw()
if (length(unique(dfPlot$Cluster)) > 15){
plotList[[tag]] <- plotList[[tag]] + theme(legend.position = "none")
}
FNbase <- paste0(tag, VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotList[[tag]])
dev.off()
figLegend <- paste0(
'**Figure ',
figureCount,
':** ',
' UMAP showing all cells from all samples together with the estimated cell-cycle phase color-coded. Download a pdf of this figure <a href="',FNrel,'" target="_blank">here</a>.'
)
figureCount <- figureCount + 1
NewChnk <- paste0(
"#### ", tag,
"\n```{r CC_UMAP_",
tag,", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",
figLegend,"'}\n",
"\n",
"\n print(plotList[['",tag,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
chnkVec <- c(
chnkVec,
NewChnk
)
## Done first umap all samples ##
###############################################################################
###############################################################################
## First tsne all samples together ##
tag <- paste0("tSNE_All_Samples")
dfPlot <- OsC@meta.data
pos <- grep("included", names(dfPlot))
if (length(pos) == 0){
dfPlot[["included"]] <- "+"
}
dfPlot[["cellID"]] <- row.names(dfPlot)
dfPlot$tSNE_1 <- NULL
dfPlot$tSNE_2 <- NULL
## Get tsNE coordinates ##
coord <- data.frame(OsC@reductions$tsne@cell.embeddings)
coord[["cellID"]] <- row.names(coord)
coord <-coord[coord$cellID %in% dfPlot$cellID, ]
dfPlot <- merge(dfPlot, coord, by.x = "cellID", by.y="cellID", all=T)
dfPlot[is.na(dfPlot)] <- 0
dfPlot <- dfPlot[dfPlot$tSNE_1 != 0 & dfPlot$tSNE_2 != 0,]
## Add cluster colors ##
#dfPlot[["Cluster"]] <- paste0("C", dfPlot$seurat_clusters)
#clusterVec <- as.vector(paste0("C", unique(sort(dfPlot$seurat_clusters))))
maxX <- 1.1*max(dfPlot$tSNE_1, na.rm = T)
minX <- 1.1*min(dfPlot$tSNE_1, na.rm = T)
maxY <- 1.1*max(dfPlot$tSNE_2, na.rm = T)
minY <- 1.1*min(dfPlot$tSNE_2, na.rm = T)
#library(scales)
#clusterCols = hue_pal()(length(clusterVec))
#dfPlot$Cluster <- factor(dfPlot$Cluster, levels = clusterVec)
# dotsize = 1.5
# if (nrow(dfPlot) > 10000){
# dotsize = 0.75
# } else if (nrow(dfPlot) > 50000){
# dotsize = 0.5
# } else {
# dotsize = 0.25
# }
plotList[[tag]] <- ggplot(data=dfPlot[dfPlot$included == "+",], aes(tSNE_1, tSNE_2, color=Phase)
) + geom_point( shape=16, size = as.numeric(dotsize)
) + xlab("tSNE1") + ylab("tSNE2") + theme(
axis.text.y = element_text(size=8),
axis.text.x = element_text(size=8),
axis.title.y = element_text(size=8),
axis.title.x = element_text(size=8),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title = element_text(hjust = 0.5, size = 12),
legend.title = element_blank()
) + ggtitle(paste0("Sample: ", tag)
) + xlim(minX, maxX) + ylim(minY, maxY
) + coord_fixed(ratio=1
) + theme_bw()
if (length(unique(dfPlot$Cluster)) > 15){
plotList[[tag]] <- plotList[[tag]] + theme(legend.position = "none")
}
FNbase <- paste0(tag, VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotList[[tag]])
dev.off()
figLegend <- paste0(
'**Figure ',
figureCount,
':** ',
' tSNE showing all cells from all samples together. The esimated cell-cycle phase is color coded. Download a pdf of this figure <a href="',FNrel,'" target="_blank">here</a>.'
)
figureCount <- figureCount + 1
NewChnk <- paste0(
"#### ", tag,
"\n```{r CC_tSNE_",
tag,", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",
figLegend,"'}\n",
"\n",
"\n print(plotList[['",tag,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
chnkVec <- c(
chnkVec,
NewChnk
)
## Done first tsne all samples ##
###############################################################################
###############################################################################
## Make one UMAP plot per sample ##
sampleVec <- sort(unique(OsC@meta.data$sampleID))
dfPlot <- OsC@meta.data
pos <- grep("included", names(dfPlot))
if (length(pos) == 0){
dfPlot[["included"]] <- "+"
}
dfPlot[["cellID"]] <- row.names(dfPlot)
## Get UMAP coordinates ##
coord <- data.frame(OsC@reductions$umap@cell.embeddings)
coord[["cellID"]] <- row.names(coord)
coord <-coord[coord$cellID %in% dfPlot$cellID, ]
dfPlot$UMAP_1 <- NULL
dfPlot$UMAP_2 <- NULL
dfPlot <- merge(dfPlot, coord, by.x = "cellID", by.y="cellID", all=T)
dfPlot[is.na(dfPlot)] <- 0
dfPlot <- dfPlot[dfPlot$UMAP_1 != 0 & dfPlot$UMAP_2 != 0,]
## Add cluster colors ##
#dfPlot[["Cluster"]] <- paste0("C", dfPlot$seurat_clusters)
#clusterVec <- as.vector(paste0("C", unique(sort(dfPlot$seurat_clusters))))
#library(scales)
#clusterCols = hue_pal()(length(clusterVec))
#dfPlot$Cluster <- factor(dfPlot$Cluster, levels = clusterVec)
maxX <- 1.1*max(dfPlot$UMAP_1, na.rm = T)
minX <- 1.1*min(dfPlot$UMAP_1, na.rm = T)
maxY <- 1.1*max(dfPlot$UMAP_2, na.rm = T)
minY <- 1.1*min(dfPlot$UMAP_2, na.rm = T)
for (i in 1:length(sampleVec)){
tag <- paste0("UMAP_CC_plot_by_", sampleVec[i])
dfPlotSel <- dfPlot[dfPlot$sampleID == sampleVec[i], ]
plotList[[tag]] <- ggplot(data=dfPlotSel[dfPlot$included == "+",], aes(UMAP_1, UMAP_2, color=Phase)
) + geom_point( shape=16, size = as.numeric(dotsize)
) + xlab("UMAP1") + ylab("UMAP2") + theme(
axis.text.y = element_text(size=8),
axis.text.x = element_text(size=8),
axis.title.y = element_text(size=8),
axis.title.x = element_text(size=8),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title = element_text(hjust = 0.5, size = 12),
legend.title = element_blank()
) + ggtitle(paste0("Sample: ", tag)
) + xlim(minX, maxX) + ylim(minY, maxY
) + coord_fixed(ratio=1
) + theme_bw()
if (length(unique(dfPlot$Cluster)) > 15){
plotList[[tag]] <- plotList[[tag]] + theme(legend.position = "none")
}
FNbase <- paste0(tag, VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotList[[tag]])
dev.off()
figLegend <- paste0(
'**Figure ',
figureCount,
':** ',
' Sample-level UMAPs. Estimated cell-cylce phase color-coded. Download a pdf of this figure <a href="',FNrel,'" target="_blank">here</a>.'
)
figureCount <- figureCount + 1
NewChnk <- paste0(
paste("#### ", tag),
"\n```{r CC_UMAP_",
tag,", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",
figLegend,"'}\n",
"\n",
"\n print(plotList[['",tag,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
chnkVec <- c(
chnkVec,
NewChnk
)
}
## Done making one umap plot per sample ##
###############################################################################
###############################################################################
## Add cluster dendrogram by sample ##
if (length(unique(OsC@meta.data$sampleID)) > 3){
library(ggtree)
Idents(OsC) <- "sampleID"
OsC <- BuildClusterTree(OsC)
tag <- paste0("Sample_Dendrogram")
OsC@tools$BuildClusterTree$tip.label <- paste0("C", OsC@tools$BuildClusterTree$tip.label)
plotList[[tag]] <- ggplot(OsC@tools$BuildClusterTree
) + geom_tree(
) + theme_tree(
) + geom_tiplab(
) + labs(title=tag
) + theme(
panel.border = element_rect(colour = "black", fill=NA, size=1),
axis.title.x=element_blank(),
plot.title = element_text(hjust = 0.5, size = 12)
) + xlim(0,1.2*max(OsC@tools$BuildClusterTree[[2]]))
## Save to file ##
FNbase <- paste0(tag,".", VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotList[[tag]])
dev.off()
figLegend <- paste0(
'**Figure ',
figureCount,
':** ',
' Clusterplot dendrogram by sample ID. ','A pdf of this figure can be downloaded <a href="',FNrel,'", target="_blank>here</a>.'
)
NewChnk <- paste0(
"#### SampleID Dendrogram",
"\n```{r ", tag, "results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",
figLegend,"'}\n",
"\n",
"\n print(plotList[['",tag,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
chnkVec <- c(
chnkVec,
NewChnk
)
figureCount <- figureCount + 1
}
## Done by sample ##
###############################################################################
###############################################################################
## Add cluster dendrogram by cluster ##
library(ggtree)
Idents(OsC) <- "seurat_clusters"
OsC <- BuildClusterTree(OsC)
tag <- paste0("Cluster_Dendrogram")
OsC@tools$BuildClusterTree$tip.label <- paste0("C", OsC@tools$BuildClusterTree$tip.label)
plotList[[tag]] <- ggplot(OsC@tools$BuildClusterTree
) + geom_tree(
) + theme_tree(
) + geom_tiplab(
) + labs(title=tag
) + theme(
panel.border = element_rect(colour = "black", fill=NA, size=1),
axis.title.x=element_blank(),
plot.title = element_text(hjust = 0.5, size = 12)
) + xlim(0,1.2*max(OsC@tools$BuildClusterTree[[2]]))
#+ xlim(-1,1.2*max(OsC@tools$BuildClusterTree$edge))
## Save to file ##
FNbase <- paste0(tag,".", VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotList[[tag]])
dev.off()
figLegend <- paste0(
'**Figure ',
figureCount,
':** ',
' Clusterplot dendrogram. ','A pdf of this figure can be downloaded <a href="',FNrel,'" target="_blank">here</a>.'
)
NewChnk <- paste0(
"#### Cluster Dendrogram",
"\n```{r ", tag, "results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",
figLegend,"'}\n",
"\n",
"\n print(plotList[['",tag,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
chnkVec <- c(
chnkVec,
NewChnk
)
figureCount <- figureCount + 1
## Done integraed analysis ##
###############################################################################cat(paste(knit(text = chnkVec, quiet = T), collapse = '\n'))CellCyclePhase_All_Samples
Figure 7: UMAP showing all cells from all samples together with the estimated cell-cycle phase color-coded. Download a pdf of this figure here.
tSNE_All_Samples
Figure 8: tSNE showing all cells from all samples together. The esimated cell-cycle phase is color coded. Download a pdf of this figure here.
UMAP_CC_plot_by_S
Figure 9: Sample-level UMAPs. Estimated cell-cylce phase color-coded. Download a pdf of this figure here.
Cluster Dendrogram
Figure 10: Clusterplot dendrogram. A pdf of this figure can be downloaded here.
Barchart N cells and Percent in Clusters
if you could prepare a bar-shape graph with the %of cells clusters representing our populations (like in the Nat Med)
###############################################################################
## Create datatable for plotting ##
## This plotting procedure requires three sets: the sets cellTypeIDs, clusterIDs, cellTypeIDs
## level1ID, level2ID, level3ID
sampleIDs <- unique(OsC@meta.data$sampleID)
#Obio@parameterList$singleCellSeuratMtCutoff <- 20
if (is.null(Obio@parameterList$singleCellSeuratMtCutoff)){
Obio@parameterList$singleCellSeuratMtCutoff <- rep(10, length(sampleIDs))
} else if (length(Obio@parameterList$singleCellSeuratMtCutoff) == 1){
Obio@parameterList$singleCellSeuratMtCutoff <- rep(
Obio@parameterList$singleCellSeuratMtCutoff,
length(sampleIDs)
)
} else if (length(Obio@parameterList$singleCellSeuratMtCutoff) == length(sampleIDs)){
Obio@parameterList$singleCellSeuratMtCutoff <- Obio@parameterList$singleCellSeuratMtCutoff
} else {
print("Can't determine mitochondrial cut off. ")
stop()
}
sampleIDs <- unique(OsC@meta.data$sampleID)
clusterIDs <- unique(OsC@meta.data[,Obio@parameterList$singleCellClusterString])
if (length(grep("cellIdent", names(OsC@meta.data))) == 0){
OsC@meta.data[["cellIdent"]] <- "All"
}
cellTypeIDs <- unique(OsC@meta.data[,"cellIdent"])
dfTemp <- OsC@meta.data
if (length(grep("^cellIdent$", names(dfTemp))) == 0){
dfTemp[["cellIdent"]] <- "All"
}
dfTemp <- dfTemp[dfTemp$percent.mt <= max(Obio@parameterList$singleCellSeuratMtCutoff), ]
dfTemp[["cellID"]] <- row.names(dfTemp)
dfTemp <- unique(dfTemp[,c("cellID", "sampleID", Obio@parameterList$singleCellClusterString, "cellIdent")])
names(dfTemp) <- gsub(Obio@parameterList$singleCellClusterString, "Cluster", names(dfTemp) )
names(dfTemp) <- gsub(Obio@parameterList$singleCellClusterString, "Cluster", names(dfTemp) )
dfTemp <- unique(dfTemp[,c("cellID", "sampleID", "Cluster", "cellIdent")])
dfRes <- dfTemp
dfRes$cellID <- NULL
row.names(dfRes) <- NULL
dfRes <- unique(dfRes)
dfRes[["N_cells"]] <- 0
for (i in 1:nrow(dfRes)){
dfRes[i, "N_cells"] <- nrow(dfTemp[dfTemp$sampleID == dfRes[i,"sampleID"] & dfTemp$Cluster == dfRes[i,"Cluster"] & dfTemp$cellIdent == dfRes[i,"cellIdent"], ])
}
## Calculate cluster percentages per celltypeID ##
dfRes[["Perc_cells"]] <- 0
for (i in 1:length(cellTypeIDs)){
dfResTemp2 <- dfRes[dfRes$cellIdent == cellTypeIDs[i], ]
tempCluster <- as.vector(unique(dfResTemp2$Cluster))
for (j in 1:length(tempCluster)){
dfResTemp3 <- dfResTemp2[dfResTemp2$Cluster == tempCluster[j],]
NclusterTotal <- sum(dfResTemp3[, "N_cells"])
dfResTemp3[,"Perc_cells"] <- round(dfResTemp3[,"N_cells"]/NclusterTotal, 4)*100
if (j ==1){
dfRes3 <- dfResTemp3
} else {
dfRes3 <- rbind(dfResTemp3, dfRes3)
}
}
if (i ==1){
dfRes4 <- dfRes3
} else {
dfRes4 <- rbind(dfRes3, dfRes4)
}
}
dfRes <- dfRes4
plotListNcells <- list()
plotListPercent <- list()
chnkVec <- as.vector(NULL, mode = "character")
for (i in 1:length(cellTypeIDs)){
#############################################################################
## Create cell number plot ##
tag <- cellTypeIDs[i]
dfResTemp <- dfRes[dfRes$cellIdent == cellTypeIDs[i], ]
## Calculate percentages for this subset ##
plotListNcells[[tag]] <- ggplot(
) + geom_bar(aes(y = N_cells, x = Cluster, fill = sampleID), data = dfResTemp, stat="identity"
) + labs(title=tag, x="Cluster", y = "Cell Count"
) + theme(
panel.border = element_rect(colour = "black", fill=NA, size=1),
axis.title.x=element_blank(),
plot.title = element_text(hjust = 0.5, size = 12)
) + coord_flip()
###########################################################################
## Save plot to file ##
FNbase <- paste0(tag,".Ncells", VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotListNcells[[tag]])
dev.off()
## ##
###########################################################################
###########################################################################
## Add to chunk ##
figCap <- paste0(
'**Figure ',
figureCount,
'A:** Cell Count in each cluster for ',
tag,
'Download a pdf of this figure <a href="',FNrel,'" target="_blank">here</a>. '
)
NewChnk <- paste0(
paste0("#### Barchart_ ", tag),
"\n```{r Barchart-",tag,", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",figCap,"'}\n",
"\n",
"\n print(plotListNcells[['",tag,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
chnkVec <- c(
chnkVec,
NewChnk
)
## Done adding ##
#############################################################################
#############################################################################
## Add percentage plot ##
plotListPercent[[tag]] <- ggplot(
) + geom_bar(aes(x = Cluster, y = Perc_cells, fill = sampleID), data = dfResTemp, stat="identity"
) + labs(title=tag, x="Cluster", y = "Percent Cells"
) + theme(
panel.border = element_rect(colour = "black", fill=NA, size=1),
axis.title.x=element_blank(),
plot.title = element_text(hjust = 0.5, size = 12)
) + coord_flip()
###########################################################################
## Save plot to file ##
FNbase <- paste0(tag, ".percent.cells",VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotListPercent[[tag]])
dev.off()
## ##
###########################################################################
###########################################################################
## Add to chunk ##
figCap <- paste0(
'**Figure ',
figureCount,
'B:** Cell percentages in each cluster for ',
tag,
'Download a pdf of this figure <a href="',FNrel,'" target="_blank">here</a>. '
)
NewChnk <- paste0(
"\n```{r Barchart-percent_",tag,", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",figCap,"'}\n",
"\n",
"\n print(plotListPercent[['",tag,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
chnkVec <- c(
chnkVec,
NewChnk
)
## Done adding percentage plot ##
#############################################################################
figureCount <- figureCount + 1
}
## Done creating data table ##
###############################################################################cat(paste(knit(text = chnkVec, quiet = T), collapse = '\n'))
Category Enrichment Scatterplots
# save(Obio,
# file = paste0(
# Obio@parameterList$localWorkDir,
# Obio@parameterList$project_id,
# ".temp.bioLOGIC.Robj"
# )
# )
#print("Obio Object saved.")
save(OsC,
file = paste0(
Obio@parameterList$localWorkDir,
Obio@parameterList$project_id,
".Seurat.Robj"
)
)
library(AUCell)
plotList <- list()
chnkVec <- as.vector(NULL, mode = "character")
# Defined in the section above #
## This needs to become a gmt file ##
if (is.null(Obio@parameterList$catRefFile)){
FNcat <- paste0(hpc.mount, "Projects/schaefera/tobias.ackels/360_scRNAseq_mm_10X_1M_neurons_20k/basedata/asl320.referenceCats.txt")
} else {
FNcat <- Obio@parameterList$catRefFile
}
if (length(grep(".gmt$", FNcat)) > 0){
print("Load gmt file. To be implemented.")
stop()
} else {
dfHeatmapGenes <- read.delim(
FNcat,
header = T,
sep = "\t",
stringsAsFactors = F
)
if (is.null(Obio@parameterList[["cat2DotplotList"]])){
Obio@parameterList[["cat2DotplotList"]] <- list()
}
if (is.null(Obio@parameterList[["cat2HMplotList"]])){
Obio@parameterList[["cat2HMplotList"]] <- list()
}
for (i in 1:ncol(dfHeatmapGenes)){
genes <- as.vector(dfHeatmapGenes[2:nrow(dfHeatmapGenes),i])
genes <- genes[genes %in% rownames(x = OsC@assays$RNA)]
if (length(unique(genes)) < 61 | (length(unique(genes)) > 0)){
Obio@parameterList[["cat2DotplotList"]][[names(dfHeatmapGenes)[i]]] <- genes
}
if ((length(unique(genes)) < 501) | (length(unique(genes)) > 2) ){
Obio@parameterList[["cat2HMplotList"]] [[names(dfHeatmapGenes)[i]]] <- genes
}
}
}
## Add transcription factors to dotplot ##
if (Obio@parameterList$geneIDcolumn != "mgi_symbol" & Obio@parameterList$geneIDcolumn != "hgnc_symbol") {
queryGS <- "hgnc_symbol"
} else {
queryGS <- Obio@parameterList$geneIDcolumn
}
tempVec <- retrieve.gene.category.from.db(
cat_id = "ag_lab_categories__10",
password = db.pwd,
gene.symbol = queryGS,
user = Obio@parameterList$db.user,
host = Obio@parameterList$host
)
###############################################################################
## If this is fish, translation is non-human or non-mouse, translation is necessary
if (queryGS != Obio@parameterList$geneIDcolumn){
dfAnno <- Obio@dfGeneAnnotation
dfAnno <- unique(dfAnno[,c("hgnc_symbol",Obio@parameterList$geneIDcolumn )])
dfAnno <- dfAnno[dfAnno$hgnc_symbol != "", ]
dfAnno <- dfAnno[dfAnno$hgnc_symbol %in% tempVec, ]
tempVec <- unique(dfAnno[,Obio@parameterList$geneIDcolumn])
tempVec <- tempVec[tempVec != ""]
}
dfHMG <- dfGeneralMarkers[dfGeneralMarkers$gene %in% tempVec, ]
dfHMGsel <- data.frame(dfHMG %>% group_by(cluster) %>% top_n(5, avg_diff))
Obio@parameterList[["cat2DotplotList"]][["Top5_TF_per_cluster_Markers"]] <- as.vector(unique(dfHMGsel$gene))
## Add cluster defining transcription factors to the collection ##
## For the dotplot ##
###############################################################################
## Get backdrop
exprMatrix <- as.matrix(OsC@assays$RNA@counts)
#logMat <- log10(exprMatrix+1)
# When using a Seurat object #
logMat <- data.frame(OsC[["RNA"]]@data)
## Load tSNE coordinates ##
cellsTsne <- data.frame(OsC@reductions$umap@cell.embeddings)
## done
FNbase <- paste0("CatScatter_Rankings", VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
cells_rankings <- AUCell_buildRankings(exprMatrix)
dev.off()
geneSets <- Obio@parameterList$cat2DotplotList
cells_AUC <- AUCell_calcAUC(geneSets, cells_rankings, aucMaxRank=nrow(cells_rankings)*0.05)
## Select thresholds ##
FNbase <- paste0("CatScatterHist", VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
set.seed(123)
cells_assignment <- AUCell_exploreThresholds(
cells_AUC,
plotHist=TRUE,
nCores=1,
assign=TRUE
)
dev.off()
## Add data to dfExpr ##
## Plot CatScatters ##
for (i in 1:length(Obio@parameterList$cat2DotplotList)){
HMname <- names(Obio@parameterList$cat2DotplotList)[i]
tag <- gsub("[.]", "_", HMname)
FNbase <- paste0("CatScatterHist_", HMname, VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
selectedThresholds <- cells_assignment[[i]]$aucThr$thresholds
if ("minimumDens" %in% rownames(selectedThresholds)) {
pThr <- selectedThresholds["minimumDens", "threshold"]
} else if ("Global_k1" %in% rownames(selectedThresholds)){
pThr <- selectedThresholds["Global_k1", "threshold"]
} else {
pThr <- selectedThresholds[1, "threshold"]
}
if (nrow(cellsTsne) > 15000){
cex = 0.25
} else if (nrow(cellsTsne) > 1000){
cex = 0.5
} else {
cex = 1
}
## Get AUC matrix ##
tSNE.df <- data.frame(cellsTsne, cell=rownames(cellsTsne))
mAUC <- getAUC(cells_AUC)[HMname,rownames(tSNE.df)]
dfAUC <- data.frame(mAUC)
dfAUC[["cellID"]] <- row.names(dfAUC)
dfAUC <- merge(dfAUC, tSNE.df, by.x = "cellID", by.y = "cell")
dfDocAUC <- unique(dfAUC[,c("cellID", "mAUC")])
dfDocAUC[["cat"]] <- paste0("Cat_",tag)
if (i == 1){
dfResAUC <- dfDocAUC
} else {
dfResAUC <- rbind(
dfResAUC,
dfDocAUC
)
}
input <- list(
"x_axis" = "UMAP1",
"y_axis" = "UMAP2",
"gene" = HMname
)
#dotsize <- cex
legendNote <- paste0(
" The following genes of this dataset are represented in this figure: ",
paste0(sort(Obio@parameterList$cat2DotplotList[[i]]), collapse = ", ")
)
plotList[[tag]] <- ggplot(data = dfAUC, aes(x=UMAP_1, y=UMAP_2, color = mAUC)
)+ geom_point( shape=16, size = dotsize
) + scale_color_gradient("AUC", low="grey", high="darkblue"
) + xlab(input$x_axis) + ylab(input$y_axis) + theme(
axis.text.y = element_text(size=8),
axis.text.x = element_text(size=8),
axis.title.y = element_text(size=8),
axis.title.x = element_text(size=8),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title = element_text(hjust = 0.5, size = 12)
)+ ggtitle(paste0("Category: ", input$gene)
) + coord_fixed(ratio = 1
) + theme_bw()
#+ theme(legend.position="none")
FNbase <- paste0("CatScatter", HMname, VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotList[[tag]])
dev.off()
## Create R markdown chunk ##
figLegend <- paste0(
'**Figure ',
figureCount,
'A:** Category Scatter showing gene category ',
HMname, '. ', legendNote,
'. Download a pdf of this figure <a href="',FNrel,'" target="_blank">here</a>. '
)
NewChnk <- paste0(
"#### Category Feature Plot ",HMname,
"\n```{r CatFeatPlot1_",
i,", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",
figLegend,"'}\n",
"\n",
"\n print(plotList[['",tag,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
chnkVec <- c(
chnkVec,
NewChnk
)
###########################################################################
## Add part B - dotplot ##
DefaultAssay(OsC) <- "RNA"
OsC@meta.data[["hmIdent2"]] <- paste0("C", OsC@meta.data[,Obio@parameterList$singleCellClusterString])
levels <- paste0(
"C",
sort(unique(OsC@meta.data[,Obio@parameterList$singleCellClusterString]))
)
OsC@meta.data$hmIdent2 <- factor(OsC@meta.data$hmIdent2, levels=levels)
Idents(OsC) <- "hmIdent2"
HMname <- paste0("Dotplot_", names(Obio@parameterList$cat2DotplotList)[i])
tag <- gsub("[.]", "_", HMname)
dpGenes <- unique(Obio@parameterList$cat2DotplotList[[i]])
legendNote <- paste0("The following genes were found in this category and the single-cell dataset: ", paste0(dpGenes, collapse=", "))
plotList[[tag]] <- DotPlotSB(
object = OsC,
features = dpGenes,
#cols = cols,
group.by = NULL,
split.by = NULL,
dot.scale = 4,
col.min = 0,
col.max = 5,
#assay = "RNA"
) + ggtitle(gsub("_", "", tag)) + coord_fixed() + coord_flip()
FNbase <- paste0(HMname, VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotList[[tag]])
dev.off()
## Create R markdown chunk ##
figLegend <- paste0(
'**Figure ',
figureCount,
'B:** Dotplot showing gene category ',
HMname, '. ', legendNote,
'. Download a pdf of this figure <a href="',FNrel,'" target="_blank">here</a>. '
)
figureCount <- figureCount + 1
NewChnk <- paste0(
"\n```{r ",tag,
", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",
figLegend,"'}\n",
"\n",
"\n print(plotList[['",tag,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
chnkVec <- c(
chnkVec,
NewChnk
)
## Done adding dotplot ##
###########################################################################
}
Obio@dataTableList[["dfResAUC"]] <- dfResAUCcat(paste(knit(text = chnkVec, quiet = T), collapse = '\n'))Category Feature Plot ENS_glia
Figure 12A: Category Scatter showing gene category ENS_glia. The following genes of this dataset are represented in this figure: foxd3, sox10, sox2. Download a pdf of this figure here.
Figure 12B: Dotplot showing gene category Dotplot_ENS_glia. The following genes were found in this category and the single-cell dataset: foxd3, sox10, sox2. Download a pdf of this figure here.
Category Feature Plot ENS_likely_glia
Figure 13A: Category Scatter showing gene category ENS_likely_glia. The following genes of this dataset are represented in this figure: foxd3, mpz, plp1b, ptprz1a, ptprz1b, sox10, sox2. Download a pdf of this figure here.
Figure 13B: Dotplot showing gene category Dotplot_ENS_likely_glia. The following genes were found in this category and the single-cell dataset: foxd3, sox10, sox2, plp1b, mpz, ptprz1a, ptprz1b. Download a pdf of this figure here.
Category Feature Plot ENS_neuron
Figure 14A: Category Scatter showing gene category ENS_neuron. The following genes of this dataset are represented in this figure: chata, elavl3, elavl4, nmu, nos1, phox2bb, ret, vip, vipb. Download a pdf of this figure here.
Figure 14B: Dotplot showing gene category Dotplot_ENS_neuron. The following genes were found in this category and the single-cell dataset: elavl3, elavl4, ret, phox2bb, nmu, chata, nos1, vip, vipb. Download a pdf of this figure here.
Category Feature Plot ENS_ChAT_neuron
Figure 15A: Category Scatter showing gene category ENS_ChAT_neuron. The following genes of this dataset are represented in this figure: chata. Download a pdf of this figure here.
Figure 15B: Dotplot showing gene category Dotplot_ENS_ChAT_neuron. The following genes were found in this category and the single-cell dataset: chata. Download a pdf of this figure here.
Figure 16A: Category Scatter showing gene category ENS_NOS_neuron__3_5. The following genes of this dataset are represented in this figure: nos1, vip, vipb. Download a pdf of this figure here.
Figure 16B: Dotplot showing gene category Dotplot_ENS_NOS_neuron__3_5. The following genes were found in this category and the single-cell dataset: nos1, vip, vipb. Download a pdf of this figure here.
Category Feature Plot Smooth_muscle_cell
Figure 17A: Category Scatter showing gene category Smooth_muscle_cell. The following genes of this dataset are represented in this figure: anxa2b, cald1a, gata6, myh11a, mylka, srfa. Download a pdf of this figure here.
Figure 17B: Dotplot showing gene category Dotplot_Smooth_muscle_cell. The following genes were found in this category and the single-cell dataset: mylka, myh11a, cald1a, srfa, gata6, anxa2b. Download a pdf of this figure here.
Category Feature Plot Interstitial_cells_of_cajal
Figure 18A: Category Scatter showing gene category Interstitial_cells_of_cajal. The following genes of this dataset are represented in this figure: ano1, kita, kitb. Download a pdf of this figure here.
Figure 18B: Dotplot showing gene category Dotplot_Interstitial_cells_of_cajal. The following genes were found in this category and the single-cell dataset: ano1, kita, kitb. Download a pdf of this figure here.
Category Feature Plot Immune_cell
Figure 19A: Category Scatter showing gene category Immune_cell. The following genes of this dataset are represented in this figure: lck, lcp1, lyz. Download a pdf of this figure here.
Figure 19B: Dotplot showing gene category Dotplot_Immune_cell. The following genes were found in this category and the single-cell dataset: lcp1, lck, lyz. Download a pdf of this figure here.
Category Feature Plot Proliferation
Figure 20A: Category Scatter showing gene category Proliferation. The following genes of this dataset are represented in this figure: top2a. Download a pdf of this figure here.
Figure 20B: Dotplot showing gene category Dotplot_Proliferation. The following genes were found in this category and the single-cell dataset: top2a. Download a pdf of this figure here.
Category Feature Plot Apoptosis
Figure 21A: Category Scatter showing gene category Apoptosis. The following genes of this dataset are represented in this figure: casp2, casp6, casp7, casp8, casp9. Download a pdf of this figure here.
Figure 21B: Dotplot showing gene category Dotplot_Apoptosis. The following genes were found in this category and the single-cell dataset: casp2, casp6, casp7, casp8, casp9. Download a pdf of this figure here.
Category Feature Plot Zebrafish_adult_RG_Top_8_genes
Figure 22A: Category Scatter showing gene category Zebrafish_adult_RG_Top_8_genes. The following genes of this dataset are represented in this figure: atp1a1b, fabp7a, glula, her4.1, mdka, si:ch211-251b21.1, zgc:165461. Download a pdf of this figure here.
Figure 22B: Dotplot showing gene category Dotplot_Zebrafish_adult_RG_Top_8_genes. The following genes were found in this category and the single-cell dataset: fabp7a, her4.1, glula, si:ch211-251b21.1, zgc:165461, atp1a1b, mdka. Download a pdf of this figure here.
Category Feature Plot Zebrafish_adult_NBN_1_Top_5_genes
Figure 23A: Category Scatter showing gene category Zebrafish_adult_NBN_1_Top_5_genes. The following genes of this dataset are represented in this figure: cd99l2, cnp, ppp1r14ba, tubb5. Download a pdf of this figure here.
Figure 23B: Dotplot showing gene category Dotplot_Zebrafish_adult_NBN_1_Top_5_genes. The following genes were found in this category and the single-cell dataset: tubb5, cd99l2, cnp, ppp1r14ba. Download a pdf of this figure here.
Category Feature Plot Zebrafish_adult_NBN_2_Top_5_genes
Figure 24A: Category Scatter showing gene category Zebrafish_adult_NBN_2_Top_5_genes. The following genes of this dataset are represented in this figure: ebf3a, msi2b, tbr1b. Download a pdf of this figure here.
Figure 24B: Dotplot showing gene category Dotplot_Zebrafish_adult_NBN_2_Top_5_genes. The following genes were found in this category and the single-cell dataset: ebf3a, tbr1b, msi2b. Download a pdf of this figure here.
Category Feature Plot Zebrafish_adult_MN_Top_5_genes
Figure 25A: Category Scatter showing gene category Zebrafish_adult_MN_Top_5_genes. The following genes of this dataset are represented in this figure: gad1b, sprn2, sult4a1, sybu, ube2ka. Download a pdf of this figure here.
Figure 25B: Dotplot showing gene category Dotplot_Zebrafish_adult_MN_Top_5_genes. The following genes were found in this category and the single-cell dataset: sult4a1, sybu, ube2ka, gad1b, sprn2. Download a pdf of this figure here.
Category Feature Plot Zebrafish_adult_OPC_Top_5_genes
Figure 26A: Category Scatter showing gene category Zebrafish_adult_OPC_Top_5_genes. The following genes of this dataset are represented in this figure: aplnra, cd82a, sema5a, traf4a. Download a pdf of this figure here.
Figure 26B: Dotplot showing gene category Dotplot_Zebrafish_adult_OPC_Top_5_genes. The following genes were found in this category and the single-cell dataset: aplnra, sema5a, traf4a, cd82a. Download a pdf of this figure here.
Category Feature Plot Zebrafish_adult_RG_Top_30_genes
Figure 27A: Category Scatter showing gene category Zebrafish_adult_RG_Top_30_genes. The following genes of this dataset are represented in this figure: atp1a1b, atp1b4, cox4i2, cx43, dla, efhd1, fabp7a, fads2, glula, her15.1, her15.2, her4.1, her4.2, her4.2, her4.4, her6, her9, lix1, mdka, ptn, s100b, s1pr1, selenop, si:ch211-251b21.1, slc1a2b, slc1a3b. Download a pdf of this figure here.
Figure 27B: Dotplot showing gene category Dotplot_Zebrafish_adult_RG_Top_30_genes. The following genes were found in this category and the single-cell dataset: fabp7a, si:ch211-251b21.1, her4.1, her4.2, atp1a1b, slc1a2b, her6, her4.4, mdka, glula, slc1a3b, s100b, cx43, her15.2, her15.1, her9, selenop, dla, atp1b4, fads2, efhd1, cox4i2, lix1, ptn, s1pr1. Download a pdf of this figure here.
Category Feature Plot Zebrafish_adult_NBN_1_Top_30_genes
Figure 28A: Category Scatter showing gene category Zebrafish_adult_NBN_1_Top_30_genes. The following genes of this dataset are represented in this figure: cct4, cd99l2, cnp, dpysl5a, elavl3, fkbp2, fscn1a, gpm6ab, hnrnpa0l, jun, kdm6bb, marcksb, mex3b, mllt11, nme2b.1, ppp1r14ba, rcan1a, rps6, si:ch211-288g17.3, si:ch73-281n10.2, tmsb, tuba1a, tuba1c, tuba2, tubb2b, tubb4b, tubb5, vat1, zfpm2a. Download a pdf of this figure here.
Figure 28B: Dotplot showing gene category Dotplot_Zebrafish_adult_NBN_1_Top_30_genes. The following genes were found in this category and the single-cell dataset: tubb5, cd99l2, cnp, ppp1r14ba, tuba2, tmsb, tuba1a, tuba1c, tubb2b, vat1, fscn1a, hnrnpa0l, mex3b, marcksb, fkbp2, kdm6bb, mllt11, gpm6ab, rcan1a, zfpm2a, dpysl5a, tubb4b, si:ch211-288g17.3, si:ch73-281n10.2, elavl3, rps6, nme2b.1, jun, cct4. Download a pdf of this figure here.
Category Feature Plot Zebrafish_adult_NBN_2_Top_30_genes
Figure 29A: Category Scatter showing gene category Zebrafish_adult_NBN_2_Top_30_genes. The following genes of this dataset are represented in this figure: atp1b1b, c1qtnf4, camk1gb, ebf3a, edil3a, elavl4, etv5a, grm8b, hsp70.1, hsp70l, kcnj19a, khdrbs2, msi2b, oxct1a, pcp4l1, ptprn2, rtn4r, snap25b, stxbp1b, tbr1b, tmeff2a, tp53i11b, uncx, uncx4.1, zfhx4. Download a pdf of this figure here.
Figure 29B: Dotplot showing gene category Dotplot_Zebrafish_adult_NBN_2_Top_30_genes. The following genes were found in this category and the single-cell dataset: ebf3a, kcnj19a, tbr1b, msi2b, rtn4r, stxbp1b, etv5a, camk1gb, tp53i11b, atp1b1b, c1qtnf4, tmeff2a, khdrbs2, pcp4l1, uncx4.1, ptprn2, edil3a, snap25b, zfhx4, elavl4, oxct1a, uncx, hsp70l, hsp70.1, grm8b. Download a pdf of this figure here.
Category Feature Plot Zebrafish_adult_MN_Top_30_genes
Figure 30A: Category Scatter showing gene category Zebrafish_adult_MN_Top_30_genes. The following genes of this dataset are represented in this figure: adcyap1b, aldoaa, atp6v0b, atp6v1h, cyfip2, dnajc5aa, fam107b, gad1b, gas7a, gnao1a, napgb, ngb, nrgna, nsfa, pam, pcsk1nl, pik3ip1, ppp3cb, rnasekb, serpini1, si:dkey-35i13.1, sprn2, sult4a1, sybu, ube2ka, upf3a. Download a pdf of this figure here.
Figure 30B: Dotplot showing gene category Dotplot_Zebrafish_adult_MN_Top_30_genes. The following genes were found in this category and the single-cell dataset: sult4a1, sybu, ube2ka, gad1b, sprn2, atp6v1h, aldoaa, pam, gas7a, nrgna, ngb, nsfa, adcyap1b, ppp3cb, serpini1, napgb, pik3ip1, atp6v0b, fam107b, gnao1a, cyfip2, pcsk1nl, si:dkey-35i13.1, upf3a, rnasekb, dnajc5aa. Download a pdf of this figure here.
Category Feature Plot Zebrafish_adult_OPC_Top_30_genes
Figure 31A: Category Scatter showing gene category Zebrafish_adult_OPC_Top_30_genes. The following genes of this dataset are represented in this figure: aplnra, aplnrb, atp1a1b, bcan, cd63, cd82a, cd9b, cxcl12a, dscamb, fabp7a, igsf9ba, lrrn1, nptx1l, olig1, olig2, rgcc, scg3, sema5a, si:busm1-57f23.1, si:ch211-132g1.3, si:ch211-137a8.4, slc1a2b, sox10, syt9b, taok3a, tmem178b, traf4a, usp3. Download a pdf of this figure here.
Figure 31B: Dotplot showing gene category Dotplot_Zebrafish_adult_OPC_Top_30_genes. The following genes were found in this category and the single-cell dataset: aplnra, sema5a, traf4a, cd82a, si:ch211-132g1.3, cd9b, olig2, aplnrb, syt9b, cxcl12a, igsf9ba, si:busm1-57f23.1, tmem178b, slc1a2b, lrrn1, atp1a1b, si:ch211-137a8.4, olig1, bcan, rgcc, dscamb, fabp7a, cd63, sox10, nptx1l, taok3a, scg3, usp3. Download a pdf of this figure here.
Category Feature Plot Zebrafish_adult_qRGC_Top_47_genes
Figure 32A: Category Scatter showing gene category Zebrafish_adult_qRGC_Top_47_genes. The following genes of this dataset are represented in this figure: acbd7, aldocb, anxa11b, atp1a1b, atp1b4, b2ml, ckbb, cx43, cxcl12a, efhd1, fabp7a, gapdhs, glula, gmnn, grm2b, hepacama, hipk2, itm2ba, lgals2a, luzp2, mfge8a, mgll, mt2, mt-atp6, mt-atp8, mt-co2, mt-co3, mt-cyb, mt-nd1, mt-nd2, psap, ptgdsb.1, ptgdsb.2, ptn, s100b, s1pr1, sdcbp2, selenop, si:ch211-251b21.1, slc1a2b, slc4a4a, slc6a11b, tegt, zgc:153704. Download a pdf of this figure here.
Figure 32B: Dotplot showing gene category Dotplot_Zebrafish_adult_qRGC_Top_47_genes. The following genes were found in this category and the single-cell dataset: fabp7a, si:ch211-251b21.1, aldocb, mt-atp6, ptn, mt-atp8, ckbb, cx43, mt-cyb, acbd7, glula, mt2, mfge8a, slc1a2b, efhd1, selenop, ptgdsb.1, mt-co2, atp1a1b, luzp2, hipk2, ptgdsb.2, gapdhs, s100b, mt-co3, slc4a4a, sdcbp2, slc6a11b, zgc:153704, tegt, psap, s1pr1, lgals2a, gmnn, b2ml, hepacama, mt-nd1, mgll, atp1b4, anxa11b, itm2ba, mt-nd2, cxcl12a, grm2b. Download a pdf of this figure here.
Category Feature Plot Zebrafish_adult_aRGC_Top_47_genes
Figure 33A: Category Scatter showing gene category Zebrafish_adult_aRGC_Top_47_genes. The following genes of this dataset are represented in this figure: actl6a, arrdc3b, ascl1a, ascl1b, boc, ccnd1, ccnd2a, cct4, ddx21, dla, h3f3a, her15.1, her4.2, her4.4, hmga1a, hmgb2b, hnrnpabb, kdm6ba, lfng, marcksb, mcm7, mibp2, nap1l1, ran, ranbp1, rpl22l1, rps23, rps6kb1b, rsl1d1, si:ch211-222l21.1, si:ch211-288g17.3, si:ch73-46j18.5, slc6a15, snrpd1, snrpd3l, sox4a, srsf2a, stmn1b, tcp1, tmsb. Download a pdf of this figure here.
Figure 33B: Dotplot showing gene category Dotplot_Zebrafish_adult_aRGC_Top_47_genes. The following genes were found in this category and the single-cell dataset: ccnd1, si:ch211-222l21.1, snrpd1, actl6a, ranbp1, ascl1b, sox4a, rsl1d1, mcm7, si:ch211-288g17.3, rps6kb1b, lfng, ascl1a, tmsb, tcp1, hmgb2b, her4.2, si:ch73-46j18.5, marcksb, hmga1a, cct4, ccnd2a, nap1l1, srsf2a, ddx21, hnrnpabb, arrdc3b, mibp2, boc, snrpd3l, her4.4, dla, kdm6ba, slc6a15, her15.1, stmn1b, rps23, ran, h3f3a, rpl22l1. Download a pdf of this figure here.
Category Feature Plot Neural_Plate_Border_Module
Figure 34A: Category Scatter showing gene category Neural_Plate_Border_Module. The following genes of this dataset are represented in this figure: dlx5a, dlx6a, gbx1, gbx2, msx1a, msx1b, msx2a, msx2b, pax3a, pax3b, pax7a, pax7b, tfap2a, tfap2b, tfap2c, tfap2d, tfap2e, zic1, zic2a, zic2b, zic3, zic4, zic5, zic6. Download a pdf of this figure here.
Figure 34B: Dotplot showing gene category Dotplot_Neural_Plate_Border_Module. The following genes were found in this category and the single-cell dataset: zic1, zic2a, zic2b, zic3, zic4, zic5, zic6, msx1a, msx1b, msx2a, msx2b, gbx1, gbx2, pax3a, pax3b, pax7a, pax7b, tfap2a, tfap2b, tfap2c, tfap2d, tfap2e, dlx5a, dlx6a. Download a pdf of this figure here.
Category Feature Plot Neural_Plate_Border_Module.1
Figure 35A: Category Scatter showing gene category Neural_Plate_Border_Module.1. The following genes of this dataset are represented in this figure: ets1, foxd3, id1, id2a, id2b, id3, id4, myb, myca, mycla, myclb, mycn, pax3a, pax3b, pax7a, pax7b, snai1a, snai1b, snai2, snai3, sox10, sox5, sox9b, tfap2b, tfap2c, tfap2d, tfap2e, twist1a, twist1b, twist2, twist3. Download a pdf of this figure here.
Figure 35B: Dotplot showing gene category Dotplot_Neural_Plate_Border_Module.1. The following genes were found in this category and the single-cell dataset: foxd3, snai1a, snai1b, snai2, snai3, twist1a, twist1b, twist2, twist3, sox5, pax3a, pax3b, pax7a, pax7b, sox9b, ets1, myca, mycla, myclb, mycn, tfap2b, tfap2c, tfap2d, tfap2e, id1, id2a, id2b, id3, id4, sox10, myb. Download a pdf of this figure here.
Category Feature Plot Neural_Plate_Border_Module.2
Figure 36A: Category Scatter showing gene category Neural_Plate_Border_Module.2. The following genes of this dataset are represented in this figure: ebf1b, foxd3, id1, id2a, id2b, id3, id4, myb, myca, mycb, mycla, myclb, mycn, pax3a, pax3b, pax7a, pax7b, rxrga, rxrgb, snai1a, snai1b, snai2, sox10, sox5, sox9a, sox9b, tfap2b, tfap2c, tfap2d, tfap2e. Download a pdf of this figure here.
Figure 36B: Dotplot showing gene category Dotplot_Neural_Plate_Border_Module.2. The following genes were found in this category and the single-cell dataset: sox10, sox5, foxd3, ebf1b, pax3a, pax3b, pax7a, pax7b, tfap2b, tfap2c, tfap2d, tfap2e, rxrga, rxrgb, myca, mycb, mycla, myclb, mycn, sox9a, sox9b, id1, id2a, id2b, id3, id4, snai1a, snai1b, snai2, myb. Download a pdf of this figure here.
Category Feature Plot Betancur_Neural_Plate_Border_Module_Targets
Figure 37A: Category Scatter showing gene category Betancur_Neural_Plate_Border_Module_Targets. The following genes of this dataset are represented in this figure: dlx5a, gbx2, her9, msx1a, msx1b, msx2a, msx2b, pax3a, pax3b, pax7a, pax7b, tfap2a, zic1, zic5. Download a pdf of this figure here.
Figure 37B: Dotplot showing gene category Dotplot_Betancur_Neural_Plate_Border_Module_Targets. The following genes were found in this category and the single-cell dataset: zic1, zic5, msx1a, msx1b, msx2a, msx2b, tfap2a, dlx5a, pax3a, pax3b, pax7a, pax7b, her9, gbx2. Download a pdf of this figure here.
Category Feature Plot Betancur_Pre.migratory_CNCC_Targets
Figure 38A: Category Scatter showing gene category Betancur_Pre.migratory_CNCC_Targets. The following genes of this dataset are represented in this figure: cdh1, cdh2, cdh6, ets1, foxd3, id1, id2a, id2b, id3, id4, myca, mycla, myclb, mycn, pax3a, pax3b, pax7a, pax7b, snai1a, snai1b, snai2, sox10, sox5, sox9a, sox9b, tfap2a, tfap2b, tfap2c, tfap2d, tfap2e, twist1a, twist1b, twist2, twist3. Download a pdf of this figure here.
Figure 38B: Dotplot showing gene category Dotplot_Betancur_Pre.migratory_CNCC_Targets. The following genes were found in this category and the single-cell dataset: foxd3, snai1a, snai1b, snai2, twist1a, twist1b, twist2, twist3, sox5, pax3a, pax3b, pax7a, pax7b, sox9a, sox9b, ets1, myca, mycla, myclb, mycn, tfap2a, tfap2b, tfap2c, tfap2d, tfap2e, id1, id2a, id2b, id3, id4, sox10, cdh1, cdh6, cdh2. Download a pdf of this figure here.
Category Feature Plot Betancur_Migratory_CNCC_Targets
Figure 39A: Category Scatter showing gene category Betancur_Migratory_CNCC_Targets. The following genes of this dataset are represented in this figure: cdh7a, cdh7b, foxd3, id1, id2a, id2b, id3, id4, itgb1a, itgb1b, nfkb1, nfkb2, pax3a, pax3b, pax7a, pax7b, snai1a, snai1b, snai2, sox10, sox9a, sox9b, tfap2a. Download a pdf of this figure here.
Figure 39B: Dotplot showing gene category Dotplot_Betancur_Migratory_CNCC_Targets. The following genes were found in this category and the single-cell dataset: itgb1a, itgb1b, cdh7a, cdh7b, pax3a, pax3b, pax7a, pax7b, tfap2a, id1, id2a, id2b, id3, id4, sox9a, sox9b, sox10, foxd3, snai1a, snai1b, snai2, nfkb1, nfkb2. Download a pdf of this figure here.
Category Feature Plot Betancur_Glia_targets
Figure 40A: Category Scatter showing gene category Betancur_Glia_targets. The following genes of this dataset are represented in this figure: cx31.7, gfap, her4.1, her4.2, her4.3, her4.4, her6, mbpa, mbpb, mitfa, mitfb, pmp22a, pmp22b. Download a pdf of this figure here.
Figure 40B: Dotplot showing gene category Dotplot_Betancur_Glia_targets. The following genes were found in this category and the single-cell dataset: her6, her4.1, her4.2, her4.3, her4.4, mitfa, mitfb, gfap, mbpa, mbpb, cx31.7, pmp22a, pmp22b. Download a pdf of this figure here.
Category Feature Plot Betancur_Neural_targets
Figure 41A: Category Scatter showing gene category Betancur_Neural_targets. The following genes of this dataset are represented in this figure: ascl1a, ascl1b, isl1, lhx1a, lhx1b, lhx2a, lhx2b, neurog1, olig1, olig2. Download a pdf of this figure here.
Figure 41B: Dotplot showing gene category Dotplot_Betancur_Neural_targets. The following genes were found in this category and the single-cell dataset: isl1, neurog1, ascl1a, ascl1b, olig1, olig2, lhx1a, lhx1b, lhx2a, lhx2b. Download a pdf of this figure here.
Category Feature Plot Zebrafish_G2_M_phase_cell_cycle_genes
Figure 42A: Category Scatter showing gene category Zebrafish_G2_M_phase_cell_cycle_genes. The following genes of this dataset are represented in this figure: anln, anp32e, aurka, aurkb, birc5a, birc5b, bub1, cbx5, ccnb2, cdc20, cdca8, cdk1, cenpe, cenpf, ckap2l, ckap5, cks1b, cks2, ctcf, dlgap5, ect2, g2e3, gas2l3, gtse1, hmgb2a, hmgb2b, hmmr, kif11, kif20ba, kif20bb, kif23, kif2c, lbr, mki67, ncapd2, ndc80, nek2, nuf2, nusap1, rangap1a, rangap1b, smc4, tacc3, tmpoa, tmpob, top2a, tpx2, ttk, tubb4b, ube2c. Download a pdf of this figure here.
Figure 42B: Dotplot showing gene category Dotplot_Zebrafish_G2_M_phase_cell_cycle_genes. The following genes were found in this category and the single-cell dataset: hmgb2a, hmgb2b, cdk1, nusap1, ube2c, birc5a, birc5b, tpx2, top2a, ndc80, cks2, nuf2, cks1b, mki67, tmpoa, tmpob, cenpf, tacc3, smc4, ccnb2, ckap2l, aurkb, bub1, kif11, anp32e, tubb4b, gtse1, kif20ba, kif20bb, cdc20, ttk, kif2c, rangap1a, rangap1b, ncapd2, dlgap5, cdca8, ect2, kif23, hmmr, aurka, anln, lbr, ckap5, cenpe, ctcf, nek2, g2e3, gas2l3, cbx5. Download a pdf of this figure here.
Category Feature Plot Zebrafish_S_phase_cell_cycle_genes
Figure 43A: Category Scatter showing gene category Zebrafish_S_phase_cell_cycle_genes. The following genes of this dataset are represented in this figure: atad2, blm, brip1, ccne2, cdc45, cdc6, cdca7a, cdca7b, chaf1b, dscc1, dtl, e2f8, exo1, fen1, gins2, gmnn, hells, mcm2, mcm4, mcm5, mcm6, msh2, nasp, pcna, pola1, pold3, prim1, rad51, rad51ap1, rfc2, rpa2, rrm1, rrm2, slbp, tipin, tyms, ubr7, uhrf1, unga, ungb, usp1, wdr76. Download a pdf of this figure here.
Figure 43B: Dotplot showing gene category Dotplot_Zebrafish_S_phase_cell_cycle_genes. The following genes were found in this category and the single-cell dataset: mcm5, pcna, tyms, fen1, mcm2, mcm4, rrm1, unga, ungb, gins2, mcm6, cdca7a, cdca7b, dtl, prim1, uhrf1, hells, rfc2, rpa2, nasp, rad51ap1, gmnn, wdr76, slbp, ccne2, ubr7, pold3, msh2, atad2, rad51, rrm2, cdc45, cdc6, exo1, tipin, dscc1, blm, usp1, pola1, chaf1b, brip1, e2f8. Download a pdf of this figure here.
Category Feature Plot Top5_TF_per_cluster_Markers
Figure 44A: Category Scatter showing gene category Top5_TF_per_cluster_Markers. The following genes of this dataset are represented in this figure: barx2, btg2, casz1, crema, dachd, egr1, elf3, etv1, fosb, fosl1a, foxo1a, foxo1b, foxo3b, foxp4, gata3, hand2, her6, hivep2a, hmgb2a, jun, klf12b, klf13, klf7a, lzts1, mef2aa, nfia, nr3c2, nr4a3, nr5a2, pbx1a, phox2a, phox2bb, pknox2, prox1a, rarab, rbpjl, roraa, satb2, smad9, sox5, tbx2a, tcf12, tcf7l2, trps1, usf1, zeb1b, zeb2b, zfhx3, zfpm1. Download a pdf of this figure here.
Figure 44B: Dotplot showing gene category Dotplot_Top5_TF_per_cluster_Markers. The following genes were found in this category and the single-cell dataset: etv1, gata3, nr3c2, zeb2b, nfia, rarab, tbx2a, tcf12, foxp4, fosb, foxo1a, fosl1a, her6, dachd, foxo1b, lzts1, klf7a, satb2, tcf7l2, casz1, smad9, zeb1b, pknox2, barx2, pbx1a, crema, hivep2a, sox5, trps1, mef2aa, klf13, klf12b, usf1, nr5a2, zfpm1, prox1a, elf3, hand2, phox2bb, phox2a, hmgb2a, jun, btg2, egr1, nr4a3, roraa, rbpjl, foxo3b, zfhx3. Download a pdf of this figure here.
Linear Regression
## Get average Expression for all cells
dfAvgExprAllCells <- Obio@dataTableList$dfAvglg10ExprAll
## Get average expression per cluster
dfAvgExprByCluster <- Obio@dataTableList$dfAvglg10ExprPerCluster
## Merge
dfDataTable <- merge(
dfAvgExprAllCells,
dfAvgExprByCluster,
by.x = "gene",
by.y = "gene",
all =T
)
dfDataTable[is.na(dfDataTable)] <- 0
clusterVec <- names(dfDataTable)
clusterVec <- clusterVec[!(clusterVec %in% c("gene", "all"))]
slopeVec <- as.vector(NULL, mode = "numeric")
intersectVec <- as.vector(NULL, mode = "numeric")
for (i in 1:length(clusterVec)){
LMformula <- as.formula(paste0(clusterVec[i]," ~ all"))
clusterFit <- lm(data=dfDataTable, formula = LMformula)
slopeVec <- c(
slopeVec,
coef(clusterFit)[[2]]
)
intersectVec <- c(
intersectVec,
coef(clusterFit)[[1]]
)
residuals <- round(clusterFit$residuals, 3)
dfTempResiduals <- data.frame(dfDataTable$gene, residuals)
names(dfTempResiduals) <- c("gene", paste0(clusterVec[i], "_Residuals"))
if (i ==1){
dfClusterResiduals <- dfTempResiduals
} else {
dfClusterResiduals <- merge(
dfClusterResiduals,
dfTempResiduals,
by.x = "gene",
by.y = "gene",
all =T
)
}
dfClusterResiduals[is.na(dfClusterResiduals)] <- 0
}
Obio@dataTableList[["dfClusterResiduals"]] <- dfClusterResiduals
## Make gene set with Residual marker genes ##
if (Obio@parameterList$geneIDcolumn != "hgnc_symbol" & Obio@parameterList$geneIDcolumn != "mgi_symbol"){
refGeneIDcolumn <- "hgnc_symbol"
dfAnno <- Obio@dfGeneAnnotation
dfAnno <- unique(dfAnno[,c("hgnc_symbol",Obio@parameterList$geneIDcolumn )])
dfAnno <- dfAnno[dfAnno[,Obio@parameterList$geneIDcolumn] %in% dfClusterResiduals[,"gene"],]
dfClusterResiduals <- merge(
dfClusterResiduals,
dfAnno,
by.x = "gene",
by.y = Obio@parameterList$geneIDcolumn
)
dfClusterResiduals$gene <- NULL
names(dfClusterResiduals) <- gsub("hgnc_symbol", "gene",names(dfClusterResiduals))
} else {
refGeneIDcolumn <- Obio@parameterList$geneIDcolumn
}
residualClusterMarkers <- list()
DefaultcutOff <- 1
clusterVec <- names(dfClusterResiduals)
clusterVec <- clusterVec[clusterVec != "gene"]
for (i in 1:length(clusterVec)){
clusterGenes <- as.vector(sort(unique(dfClusterResiduals[dfClusterResiduals[, clusterVec[i]] > DefaultcutOff,"gene"])))
if (length(clusterGenes) > 1){
residualClusterMarkers[[paste0(clusterVec[[i]], "_Linear_Regression_Markers_", DefaultcutOff)]] <- c(
paste0(clusterVec[[i]], "_Linear_Regression_Markers_", DefaultcutOff, " Experiment ", Obio@parameterList$project_id),
clusterGenes
)
}
}
## Upload Marker Genes for this Project ##
#######################################################################
## Upload/update category by category ##
updatedCatIDs <- as.vector(NULL, mode = "character")
updatedCatNames <- as.vector(NULL, mode = "character")
for (i in 1:length(residualClusterMarkers)){
cat.name <- names(residualClusterMarkers)[i]
cat_type <- paste0("temp_cluster_marker_", Obio@parameterList$project_id)
cat.description.text <- as.vector(residualClusterMarkers[[i]][1])
gene.vec <- as.vector(
residualClusterMarkers[[i]]
)[2:length(residualClusterMarkers[[i]])]
gene.vec <- gene.vec[gene.vec != ""]
gene.vec <- sort(na.omit(gene.vec))
## Determine if cat exists ##
catID <- add.category.to.lab.reference.table.hs(
host = Obio@dbDetailList$host,
pwd = db.pwd,
user = Obio@dbDetailList$db.user,
cat.ref.db = Obio@dbDetailList$ref.cat.db,
cat.ref.db.table = Obio@parameterList$lab.categories.table,
gene.vector = gene.vec,
gene.id = refGeneIDcolumn, #options hgnc_symbol, mgi_symbol
mm.hs.conversion.file = paste0(hpc.mount, "Projects/reference_data/20160303.homologene.data.txt"),
cat_name = cat.name,
cat_type = cat_type,
data_source = paste0(Obio@parameterList$labname, " Lab"),
comments_1 = "",
comments_2 = "",
new.lab.category.table = FALSE,
cat.description.db = "internal_categories",
cat.description.db.table = "category_description",
cat.description.text = cat.description.text,
lab.name = Obio@parameterList$labname,
replaceExistingCatName = TRUE
)
updatedCatIDs <- c(
updatedCatIDs,
catID
)
updatedCatNames <- c(
updatedCatNames,
cat.name
)
} ## End dfDat loop
dfCatNameLR <- data.frame(catID = updatedCatIDs, cat_name= updatedCatNames)
###############################################################################
## Creaate Plot array ##
plotList <- list()
chnkVec <- as.vector(NULL, mode = "character")
# x-axis: all
# y-axis: intensity cluster X
# highlight: Marker genes
dfClusterResiduals <- Obio@dataTableList$dfClusterResiduals
clusterVec <- names(dfDataTable)
clusterVec <- clusterVec[!(clusterVec %in% c("gene", "all"))]
## Determine cluster colors ##
dfClustCol <- data.frame(cluster=clusterVec, clusterCol=clusterVec)
dfClustCol$clusterCol <- as.numeric(gsub("C_", "", dfClustCol$clusterCol))
dfClustCol <- dfClustCol[order(dfClustCol$clusterCol, decreasing = F),]
library(scales)
clusterCols = hue_pal()(nrow(dfClustCol))
dfClustCol[["clusterCol"]] <- clusterCols
clusterVec <- as.vector(dfClustCol$cluster)
clusterCol <- as.vector(dfClustCol$clusterCol)
for (i in 1:length(clusterVec)){
## Determine top10 markers ##
clustCol <- names(dfClusterResiduals)[grep(paste0(clusterVec[i],"_"), names(dfClusterResiduals))]
dfClusterResiduals <- dfClusterResiduals[order(dfClusterResiduals[,clustCol], decreasing = T), ]
markerVec <- as.vector(dfClusterResiduals[1:10,"gene"])
dfTempResiduals <- dfClusterResiduals[dfClusterResiduals[,clustCol] > DefaultcutOff, ]
highlightGenes <- as.vector(dfTempResiduals[,"gene"])
dfPlot <- dfDataTable[,c("gene", "all", clusterVec[i])]
dfPlot[["label"]] <- ""
dfPlot[dfPlot$gene %in% markerVec, "label"] <- dfPlot[dfPlot$gene %in% markerVec, "gene"]
## Add highlight ##
dfPlot[["ClusterMarker"]] <- ""
dfPlot[dfPlot$gene %in% highlightGenes, "ClusterMarker"] <- "+"
dfPlot[["x"]] <- dfPlot$all
dfPlot[["y"]] <- dfPlot[, clusterVec[i]]
tag <- paste0("LinearRegressionMarkers_Cluster", clusterVec[i])
## Create fitted line ##
lineFit <- lm(data=dfPlot, formula = y ~ x)
library(ggrepel)
plotList[[tag]] <- ggplot(
data = dfPlot,
aes(
x=x,
y=y, label = label, color = ClusterMarker)
#) + geom_abline(intercept = intersectVec[i], slope = slopeVec[i], linetype="dashed", color = "grey"
) + geom_smooth(method='lm', formula= y ~ x, linetype="dashed", se = T, colour = "grey"
) + geom_point( shape=16, size = dotsize
) + xlab("Average Expression All Cells") + ylab(paste0("Average Expression ", clusterVec[[i]])) + theme(
axis.text.y = element_text(size=8),
axis.text.x = element_text(size=8),
axis.title.y = element_text(size=8),
axis.title.x = element_text(size=8),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title = element_text(hjust = 0.5, size = 12)
) + ggtitle(paste0("Cluster Markers ", clusterVec[[i]])
) + scale_color_manual(values=c("#000000", clusterCol[i])
) + geom_text_repel(
) + theme_bw()
# + xlim(0, xmax) + ylim(0, ymax)
###########################################################################
## Save plot to file ##
FNbase <- paste0("Linear.regression.cluster.markers.",clusterVec[i], VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotList[[tag]])
dev.off()
## ##
###########################################################################
pos <- grep(paste0(clusterVec[i],"_"), dfCatNameLR$cat_name)
LRcatID <- ""
if (length(pos) == 1){
LRcatID <- as.vector(dfCatNameLR[pos, "catID"])
}
pos <- grep(paste0(gsub("C_", "Cluster_",clusterVec[i]),"$"), dfFAMplotIDs$cat_name)
DGEcatID <- ""
if (length(pos) == 1){
DGEcatID <- as.vector(dfFAMplotIDs[pos, "cat_id"])
}
link <- paste0('<a href="https://',urlString,'/',Obio@parameterList$project_id,'/scatterplot?x_axis=add_counts_Avg_log10_Expr_all&y_axis=norm_counts_Avg_log10_Expr_',clusterVec[i], '&cat_id=', LRcatID,'" target="_blank">here</a>')
linkDGE <- paste0('<a href="https://',urlString,'/',Obio@parameterList$project_id,'/scatterplot?x_axis=add_counts_Avg_log10_Expr_all&y_axis=norm_counts_Avg_log10_Expr_',clusterVec[i], '&cat_id=', DGEcatID,'" target="_blank">here</a>')
heatmaplink <- paste0('In order to identify cluster specific marker genes a <a href="https://biologic.crick.ac.uk/vpl362A/category-view/', LRcatID,'" target="_blank">cluster-maker gene heatmap for the linear-regression derrived marker gene set</a> might be interesting.')
heatmaplinkDGE <- paste0('In order to identify cluster specific marker genes a <a href="https://biologic.crick.ac.uk/vpl362A/category-view/', DGEcatID,'" target="_blank">cluster-maker gene heatmap for the differential-gene expression derrived marker gene set</a> might be interesting.')
figCap <- paste0(
'**Figure ',
figureCount,
':** Average gene expression intensity of all cells together versus averaged expression intensities for cluster ',clusterVec[i],'. ',
'Download a pdf of this figure <a href="',FNrel,'" target="_blank">here</a>. ',
'An interactive version of this scatterplot can be found ', link, ' with the linear regression marker genes highlighted, and ',linkDGE,' with the differential gene expression marker genes highlighted. ',
heatmaplink,
heatmaplinkDGE
)
NewChnk <- paste0(
"#### ",tag,
"\n```{r LR_marker_",tag,", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",figCap,"'}\n",
"\n",
"\n print(plotList[['",tag,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
## Histogram Part C done ##
###########################################################################
chnkVec <- c(
chnkVec,
NewChnk
)
}
## Done creating plot array ##
###############################################################################Cluster Markers By Linear Regression {.tabset .tabset-fade .tabset-pills}
## plot list will be integrated in full figure ##
cat(paste(knit(text = chnkVec, quiet = T), collapse = '\n'))LinearRegressionMarkers_ClusterC_0
Figure 45: Average gene expression intensity of all cells together versus averaged expression intensities for cluster C_0. Download a pdf of this figure here. An interactive version of this scatterplot can be found here with the linear regression marker genes highlighted, and here with the differential gene expression marker genes highlighted. In order to identify cluster specific marker genes a cluster-maker gene heatmap for the linear-regression derrived marker gene set might be interesting.In order to identify cluster specific marker genes a cluster-maker gene heatmap for the differential-gene expression derrived marker gene set might be interesting.
LinearRegressionMarkers_ClusterC_1
Figure 45: Average gene expression intensity of all cells together versus averaged expression intensities for cluster C_1. Download a pdf of this figure here. An interactive version of this scatterplot can be found here with the linear regression marker genes highlighted, and here with the differential gene expression marker genes highlighted. In order to identify cluster specific marker genes a cluster-maker gene heatmap for the linear-regression derrived marker gene set might be interesting.In order to identify cluster specific marker genes a cluster-maker gene heatmap for the differential-gene expression derrived marker gene set might be interesting.
LinearRegressionMarkers_ClusterC_2
Figure 45: Average gene expression intensity of all cells together versus averaged expression intensities for cluster C_2. Download a pdf of this figure here. An interactive version of this scatterplot can be found here with the linear regression marker genes highlighted, and here with the differential gene expression marker genes highlighted. In order to identify cluster specific marker genes a cluster-maker gene heatmap for the linear-regression derrived marker gene set might be interesting.In order to identify cluster specific marker genes a cluster-maker gene heatmap for the differential-gene expression derrived marker gene set might be interesting.
LinearRegressionMarkers_ClusterC_3
Figure 45: Average gene expression intensity of all cells together versus averaged expression intensities for cluster C_3. Download a pdf of this figure here. An interactive version of this scatterplot can be found here with the linear regression marker genes highlighted, and here with the differential gene expression marker genes highlighted. In order to identify cluster specific marker genes a cluster-maker gene heatmap for the linear-regression derrived marker gene set might be interesting.In order to identify cluster specific marker genes a cluster-maker gene heatmap for the differential-gene expression derrived marker gene set might be interesting.
LinearRegressionMarkers_ClusterC_4
Figure 45: Average gene expression intensity of all cells together versus averaged expression intensities for cluster C_4. Download a pdf of this figure here. An interactive version of this scatterplot can be found here with the linear regression marker genes highlighted, and here with the differential gene expression marker genes highlighted. In order to identify cluster specific marker genes a cluster-maker gene heatmap for the linear-regression derrived marker gene set might be interesting.In order to identify cluster specific marker genes a cluster-maker gene heatmap for the differential-gene expression derrived marker gene set might be interesting.
LinearRegressionMarkers_ClusterC_5
Figure 45: Average gene expression intensity of all cells together versus averaged expression intensities for cluster C_5. Download a pdf of this figure here. An interactive version of this scatterplot can be found here with the linear regression marker genes highlighted, and here with the differential gene expression marker genes highlighted. In order to identify cluster specific marker genes a cluster-maker gene heatmap for the linear-regression derrived marker gene set might be interesting.In order to identify cluster specific marker genes a cluster-maker gene heatmap for the differential-gene expression derrived marker gene set might be interesting.
LinearRegressionMarkers_ClusterC_6
Figure 45: Average gene expression intensity of all cells together versus averaged expression intensities for cluster C_6. Download a pdf of this figure here. An interactive version of this scatterplot can be found here with the linear regression marker genes highlighted, and here with the differential gene expression marker genes highlighted. In order to identify cluster specific marker genes a cluster-maker gene heatmap for the linear-regression derrived marker gene set might be interesting.In order to identify cluster specific marker genes a cluster-maker gene heatmap for the differential-gene expression derrived marker gene set might be interesting.
LinearRegressionMarkers_ClusterC_7
Figure 45: Average gene expression intensity of all cells together versus averaged expression intensities for cluster C_7. Download a pdf of this figure here. An interactive version of this scatterplot can be found here with the linear regression marker genes highlighted, and here with the differential gene expression marker genes highlighted. In order to identify cluster specific marker genes a cluster-maker gene heatmap for the linear-regression derrived marker gene set might be interesting.In order to identify cluster specific marker genes a cluster-maker gene heatmap for the differential-gene expression derrived marker gene set might be interesting.
LinearRegressionMarkers_ClusterC_8
Figure 45: Average gene expression intensity of all cells together versus averaged expression intensities for cluster C_8. Download a pdf of this figure here. An interactive version of this scatterplot can be found here with the linear regression marker genes highlighted, and here with the differential gene expression marker genes highlighted. In order to identify cluster specific marker genes a cluster-maker gene heatmap for the linear-regression derrived marker gene set might be interesting.In order to identify cluster specific marker genes a cluster-maker gene heatmap for the differential-gene expression derrived marker gene set might be interesting.
LinearRegressionMarkers_ClusterC_9
Figure 45: Average gene expression intensity of all cells together versus averaged expression intensities for cluster C_9. Download a pdf of this figure here. An interactive version of this scatterplot can be found here with the linear regression marker genes highlighted, and here with the differential gene expression marker genes highlighted. In order to identify cluster specific marker genes a cluster-maker gene heatmap for the linear-regression derrived marker gene set might be interesting.In order to identify cluster specific marker genes a cluster-maker gene heatmap for the differential-gene expression derrived marker gene set might be interesting.
LinearRegressionMarkers_ClusterC_10
Figure 45: Average gene expression intensity of all cells together versus averaged expression intensities for cluster C_10. Download a pdf of this figure here. An interactive version of this scatterplot can be found here with the linear regression marker genes highlighted, and here with the differential gene expression marker genes highlighted. In order to identify cluster specific marker genes a cluster-maker gene heatmap for the linear-regression derrived marker gene set might be interesting.In order to identify cluster specific marker genes a cluster-maker gene heatmap for the differential-gene expression derrived marker gene set might be interesting.
LinearRegressionMarkers_ClusterC_11
Figure 45: Average gene expression intensity of all cells together versus averaged expression intensities for cluster C_11. Download a pdf of this figure here. An interactive version of this scatterplot can be found here with the linear regression marker genes highlighted, and here with the differential gene expression marker genes highlighted. In order to identify cluster specific marker genes a cluster-maker gene heatmap for the linear-regression derrived marker gene set might be interesting.In order to identify cluster specific marker genes a cluster-maker gene heatmap for the differential-gene expression derrived marker gene set might be interesting.
LinearRegressionMarkers_ClusterC_12
Figure 45: Average gene expression intensity of all cells together versus averaged expression intensities for cluster C_12. Download a pdf of this figure here. An interactive version of this scatterplot can be found here with the linear regression marker genes highlighted, and here with the differential gene expression marker genes highlighted. In order to identify cluster specific marker genes a cluster-maker gene heatmap for the linear-regression derrived marker gene set might be interesting.In order to identify cluster specific marker genes a cluster-maker gene heatmap for the differential-gene expression derrived marker gene set might be interesting.
Heatmap and Dotplot for the top5 cluster defining genes
In this section we will try to determine what are good markers for the clusters in this experiment.
plotList <- list()
chnkVec <- as.vector(NULL, mode = "character")
###############################################################################
## Make Heatmap ##
OsC@meta.data[["hmIdent2"]] <- paste0("C", OsC@meta.data[,Obio@parameterList$singleCellClusterString])
Idents(OsC) <- "hmIdent2"
levels <- paste0(
"C",
sort(unique(OsC@meta.data[,Obio@parameterList$singleCellClusterString]))
)
levels(OsC) <- levels
## Deal with more than 5000 cells ##
if (nrow(OsC@meta.data) > 5000){
set.seed(127)
n.cells <- 5000
OsC_HM <- OsC
OsC_HM@meta.data[["HM_sel"]] <- 0
selPos <- sample(x = nrow(OsC_HM@meta.data), size = n.cells, replace = FALSE, prob = NULL)
OsC_HM@meta.data[selPos, "HM_sel"] <- 1
OsC_HM <- subset(x = OsC_HM, subset = HM_sel == 1)
subsetString <- paste0("For this heatmap 5000 cells were randomly selected from ", nrow(OsC@meta.data)," cells in the experiment. ")
} else {
OsC_HM <- OsC
subsetString <- ""
}
## Scale Data ##
allGenes <- rownames(x = OsC_HM@assays$RNA)
OsC_HM <- ScaleData(OsC_HM, verbose = FALSE, features=allGenes)
## For the moment: resetting heatmap list to keep it short ##
Obio@parameterList[["cat2HMplotList"]] <- list()
Obio@parameterList[["cat2HMplotList"]][["Top5_Cluster_Markers"]] <- as.vector(dfTop5$gene)
## Add top transcription factors for each cluster ##
## Get transcription factor genes ##
if (Obio@parameterList$geneIDcolumn != "mgi_symbol" & Obio@parameterList$geneIDcolumn != "hgnc_symbol") {
queryGS <- "hgnc_symbol"
} else {
queryGS <- Obio@parameterList$geneIDcolumn
}
tempVec <- retrieve.gene.category.from.db(
cat_id = "ag_lab_categories__10",
password = db.pwd,
gene.symbol = queryGS,
user = Obio@parameterList$db.user,
host = Obio@parameterList$host
)
###############################################################################
## If this is fish, translation is non-human or non-mouse, translation is necessary
if (queryGS != Obio@parameterList$geneIDcolumn){
dfAnno <- Obio@dfGeneAnnotation
dfAnno <- unique(dfAnno[,c("hgnc_symbol",Obio@parameterList$geneIDcolumn )])
dfAnno <- dfAnno[dfAnno$hgnc_symbol != "", ]
dfAnno <- dfAnno[dfAnno$hgnc_symbol %in% tempVec, ]
tempVec <- unique(dfAnno[,Obio@parameterList$geneIDcolumn])
tempVec <- tempVec[tempVec != ""]
}
dfHMG <- dfGeneralMarkers[dfGeneralMarkers$gene %in% tempVec, ]
dfHMGsel <- data.frame(dfHMG %>% group_by(cluster) %>% top_n(5, avg_diff))
Obio@parameterList[["cat2HMplotList"]][["Top5_TF_Cluster_Markers"]] <- as.vector(unique(dfHMGsel$gene))
## Done with translation
for (i in 1:length(Obio@parameterList[["cat2HMplotList"]])){
tag <- paste0("HMM_", names(Obio@parameterList$cat2HMplotList)[i])
textSize <- 5
plotList[[tag]] <- DoHeatmap(
object = OsC_HM,
features = Obio@parameterList[["cat2HMplotList"]][[i]],
#group.by = "hmIdent",
draw.lines =T,
label = T,
group.bar = TRUE,
slot = "scale.data",
lines.width = 2 #With of separation lines in 'cells'
#slim.col.label = TRUE,
#remove.key = removeKey
# ) + theme(legend.position = "none"
) + theme(text = element_text(size=textSize)
) + scale_fill_gradientn(colors = c("blue", "white", "red"))
## Save to file ##
FNbase <- paste0(tag, VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotList[[tag]])
dev.off()
## Create R markdown chunk ##
figLegend <- paste0(
'**Figure ',
figureCount,
':** Heatmap showing the most distinct marker genes in each cluster. ' , subsetString,
'Download a pdf of this figure <a href="',FNrel,'" target="_blank">here</a>. '
)
figureCount <- figureCount + 1
NewChnk <- paste0(
"#### Heatmap_var_genes",
"\n```{r Heatmap_", tag,
", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",
figLegend,"'}\n",
"\n",
"\n print(plotList[['",tag,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
chnkVec <- c(
chnkVec,
NewChnk
)
}
rm(OsC_HM)
## Done making Heatmap ##
###############################################################################
###############################################################################
## Make dotplot
Idents(OsC) <- "hmIdent2"
levels <- paste0(
"C",
sort(unique(OsC@meta.data[,Obio@parameterList$singleCellClusterString]))
)
levels(OsC) <- levels
dpGenes <- as.vector(unique(dfTop5$gene))
if (length(dpGenes) >= 50){
dpGenes <- as.vector(unique(dfTop1$gene))
}
#dpGenes <- rev(dpGenes[!(duplicated(dpGenes))])
tag <- paste0("Dotplot_", "Var_Genes")
textSize <- 5
plotList[[tag]] <- DotPlotSB(
object = OsC,
features = dpGenes,
#cols = cols,
group.by = NULL,
split.by = NULL,
dot.scale = 4,
col.min = 0,
col.max = 5
) + ggtitle(gsub("_", "", tag)) + coord_fixed() + coord_flip()
## Save to file ##
FNbase <- paste0(tag, VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotList[[tag]])
dev.off()
## Create R markdown chunk ##
figLegend <- paste0(
'**Figure ',
figureCount,
':** Dotplot showing showing selected marker genes. ',
'Download a pdf of this figure <a href="',FNrel,'" target="_blank">here</a>. '
)
figureCount <- figureCount + 1
NewChnk <- paste0(
"#### Dotplot Markers",
"\n```{r Dotplot_var_",
", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",
figLegend,"'}\n",
"\n",
"\n print(plotList[['",tag,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
chnkVec <- c(
chnkVec,
NewChnk
)
## Done making dotplot ##
###############################################################################
############################
## Make cat feature plot
## done making cat feature plot
################################cat(paste(knit(text = chnkVec, quiet = T), collapse = '\n'))Heatmap_var_genes
Figure 45: Heatmap showing the most distinct marker genes in each cluster. Download a pdf of this figure here.
Heatmap_var_genes
Figure 46: Heatmap showing the most distinct marker genes in each cluster. Download a pdf of this figure here.
Dotplot Markers
Figure 47: Dotplot showing showing selected marker genes. Download a pdf of this figure here.
Cluster-defining Genes Table
library(DT)
###############################################################################
## Calculate percentages of expressed genes ##
DefaultAssay(OsC) <- "RNA"
my_genes <- rownames(x = OsC@assays$RNA)
exp <- FetchData(OsC, my_genes)
ExprMatrix <- round(as.matrix(colMeans(exp > 0)) *100,1)
colnames(ExprMatrix)[1] <- "count_cut_off"
dfExprMatrix <- data.frame(ExprMatrix)
dfExprMatrix[["gene"]] <- row.names(dfExprMatrix)
Obio@dataTableList[["dfPercCellsExpr"]] <- dfExprMatrix
## Done calculating percentages of expressed gens ##
###############################################################################
###############################################################################
## Create one table per cluster ##
## Add expressed in N percent cells ##
dfPercCellsExpr <- Obio@dataTableList$dfPercCellsExpr
dfDat <- Obio@dataTableList$dfGeneralMarkersFilt
dfDat$avg_diff <- round(dfDat$avg_diff,2)
dfDat <- dfDat[,c("cluster", "gene", "avg_diff", "power")]
names(dfDat) <- gsub("avg_diff", "DGE_avg_diff", names(dfDat))
names(dfDat) <- gsub("power", "DGE_power", names(dfDat))
dfDat[["join_col"]] <- paste0(dfDat$gene, "_C_", dfDat$cluster)
dfDat$cluster <- NULL
dfDat$gene <- NULL
## Add residual results ##
dfClusterResiduals <- Obio@dataTableList$dfClusterResiduals
library(tidyverse)
dfLongResiduals <- data.frame(
dfClusterResiduals %>% pivot_longer(!gene, names_to = "cluster", values_to = "residuals")
)
###############################################################################
## Calculate Coefficient of variation for no-cluster X genes to find best unique cluster markers
## Done ##
###############################################################################
dfLongResiduals <- dfLongResiduals[abs(dfLongResiduals$residuals) > 1, ]
dfLongResiduals$cluster <- gsub("_Residuals", "", dfLongResiduals$cluster)
dfLongResiduals[["join_col"]] <- paste0(
dfLongResiduals$gene,
"_",
dfLongResiduals$cluster
)
dfLongResiduals$gene <- NULL
dfLongResiduals$cluster <- NULL
dfDat <- merge(
dfDat,
dfLongResiduals,
by.x = "join_col",
by.y = "join_col",
all =TRUE
)
dfDat[is.na(dfDat)] <- 0
names(dfDat) <- gsub("residuals", "LinearReg_Residuals", names(dfDat))
dfDat <- dfDat[order(dfDat$LinearReg_Residuals, decreasing = T),]
dfDat[["gene"]] <- sapply(
dfDat$join_col, function(x) unlist(strsplit(x, "_"))[1]
)
dfDat[["Cluster"]] <- sapply(
dfDat$join_col, function(x) unlist(strsplit(x, "_C_"))[2]
)
dfDat$join_col <- NULL
dfDat$Cluster <- paste0("Cluster_",dfDat$Cluster, "_C")
dfDat$gene <- substr(dfDat$gene,1,50)
dfDat[["uniqueMarker"]] <- as.character(!duplicated(dfDat$gene))
dfDat$uniqueMarker <- substr(dfDat$uniqueMarker, 1,1)
dtList <- list()
tabClusters <- sort(unique(dfDat$cluster))
chnkVec <- as.vector(NULL, mode="character")
linkGeneView <- paste0('<a href="https://',urlString,'/',Obio@parameterList$project_id,'/gene-view" target = "_blank">GeneView</a>')
linkFeatureView <- paste0('<a href="https://',urlString,'/mdata/',Obio@parameterList$project_id,'/html/FeatureView.html" target="_blank">FeatureView</a>')
#for (i in 1:length(tabClusters)){
#tabLegend = paste0("**Table: ** Positive and negative marker genes for ", tabClusters[i])
tabLegend = paste0("**Table: ** Positive and negative cluster-defining marker genes. Perc_Cells_Expr: Percentage of total cells expressing gene X. Enr in Cluster: Enrichment of gene X in cluster Y. To collapse the table to one particular cluster, type the name of the cluster in the search box.",
"Use the ",linkGeneView," or ",linkFeatureView," functionalities to examine individual genes in more detail. "
)
#dfTempDat <- dfDat[dfDat$cluster == tabClusters[i],]
dfTempDat <- dfDat
## Percent expressed genes
dfTempDat <- merge(
dfTempDat,
Obio@dataTableList$dfPercCellsExpr,
by.x = "gene",
by.y = "gene"
)
names(dfTempDat) <- gsub("count_cut_off", "Perc_Cells_Expr",names(dfTempDat))
names(dfTempDat) <- gsub("myAUC", "AUC", names(dfTempDat))
names(dfTempDat) <- gsub("[.]", "", names(dfTempDat))
#dtList[[paste0("Table",i)]] <- datatable(dfDat,rownames = FALSE)
if (Obio@parameterList$host == "10.27.241.234"){
urlString <- "biologic.thecrick.org"
} else {
urlString <- "biologic.crick.ac.uk"
}
dfTempDat$gene <- paste0('<a href="https://',urlString,'/',Obio@parameterList$project_id,'/gene-view?query=',dfTempDat$gene,'&exact=TRUE" target="_blank">', dfTempDat$gene, '</a>')
NewChnk <- paste0(
"#### ", names(dtList),
"\n```{r datatable_",
i,", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",
tabLegend,"'}\n",
"\n",
"\n datatable(dfTempDat,rownames = FALSE, escape = FALSE)",
"\n cat( '\n')",
"\n\n\n```\n"
)
chnkVec <- c(
chnkVec,
NewChnk
)
#}
## Done creating one table per cluster ##
##############################################################################cat(paste(knit(text = chnkVec, quiet = T), collapse = '\n'))Table: Positive and negative cluster-defining marker genes. Perc_Cells_Expr: Percentage of total cells expressing gene X. Enr in Cluster: Enrichment of gene X in cluster Y. To collapse the table to one particular cluster, type the name of the cluster in the search box.Use the GeneView or FeatureView functionalities to examine individual genes in more detail.
Define categories enriched in individual clusters
# library(knitr)
# library(ggplot2)
#
# #save.image("temp.RData")
#
# DefaultAssay(OsC) <- "RNA"
#
# dfGeneralMarkersFilt <- data.frame(
# Obio@dataTableList$dfGeneralMarkersFilt,
# stringsAsFactors = FALSE
# )
#
# dfGeneralMarkersFilt$cluster <- paste0("C", dfGeneralMarkersFilt$cluster)
#
# clusterIDs <- as.character(unique(dfGeneralMarkersFilt$cluster))
#
# plotList <- list()
# chnkVec <- as.vector(NULL, mode = "character")
#
# for (j in 1:length(clusterIDs)){
# dfEnrichSel <- dfGeneralMarkersFilt[dfGeneralMarkersFilt$cluster == as.vector(clusterIDs[j]), ]
# #dfEnrichSel[["gene"]] <- row.names(dfEnrichSel)
#
# posTestGeneSet <- as.vector(
# unique(
# dfEnrichSel[dfEnrichSel$avg_diff > lgFCthreshold, "gene"]
# )
# )
#
#
# negTestGeneSet <- as.vector(
# unique(
# dfEnrichSel[dfEnrichSel$avg_diff < -lgFCthreshold, "gene"]
# )
# )
#
# ## Get background gene set ##
# #backgroundGeneVec <- row.names(OsC[["RNA"]]@counts)
# if ((length(posTestGeneSet) >= 3) |(length(negTestGeneSet) >= 3)){
# library(enrichR)
# topMaxCat <- 10
# dbs <- listEnrichrDbs()
#
# dbs <- c("GO_Biological_Process_2017")
#
#
# PosEnriched <- enrichr(posTestGeneSet, dbs)
#
# for (i in 1:length(dbs)){
# dfTemp <- PosEnriched[[dbs[i]]]
#
# if (i ==1){
# dfPosEnriched <- dfTemp
# } else {
# dfPosEnriched <- rbind(
# dfPosEnriched,
# dfTemp
# )
# }
#
# }
#
# dfPosEnriched[["log10FDR"]] <- -1*log10(dfPosEnriched$Adjusted.P.value)
# dfPosEnriched <- dfPosEnriched[order(-dfPosEnriched$log10FDR),]
# dfPosEnriched <- na.omit(dfPosEnriched)
#
# ## Negative Side ##
# NegEnriched <- enrichr(negTestGeneSet, dbs)
#
# for (i in 1:length(dbs)){
# dfTemp <- NegEnriched[[dbs[i]]]
#
# if (i ==1){
# dfNegEnriched <- dfTemp
# } else {
# dfNegEnriched <- rbind(
# dfNegEnriched,
# dfTemp
# )
# }
#
# }
#
#
# dfNegEnriched[["log10FDR"]] <- -1*log10(dfNegEnriched$Adjusted.P.value)
# dfNegEnriched <- dfNegEnriched[order(-dfPosEnriched$log10FDR),]
# dfNegEnriched <- na.omit(dfNegEnriched)
#
# dfNegSel <- dfNegEnriched
# if (nrow(dfNegSel) > topMaxCat){
# dfNegSel <- dfNegSel[1:topMaxCat,]
# }
#
# dfPosSel <- dfPosEnriched
# if (nrow(dfPosSel) > topMaxCat){
# dfPosSel <- dfPosSel[1:topMaxCat,]
# }
#
# if ((nrow(dfNegEnriched) > 0) | (nrow(dfPosEnriched) > 0)){
#
#
# dfNegSel$log10FDR <- -1* dfNegSel$log10FDR
#
# dfSel <- rbind(
# dfNegSel,
# dfPosSel
# )
#
# dfSel <- na.omit(dfSel)
# dfSel <- dfSel[order(dfSel$log10FDR),]
# dfSel$log10FDR <- round(dfSel$log10FDR, 2)
#
# dfSel[["Category"]] <- ""
# dfSel[dfSel$log10FDR >= 0, "Category"] <- "Enr."
# dfSel[dfSel$log10FDR < 0, "Category"] <- "Depl."
#
# for (k in 1:nrow(dfSel)){
# if (nchar(dfSel[k, "Term"]) > 50 & length(grep("\\(GO", as.vector(dfSel[k, "Term"]))) > 0){
# part1 <- unlist(strsplit(as.vector(dfSel[k, "Term"]), "\\(GO"))[1]
# part1 <- substr(part1, 1, 45)
# part2 <- unlist(strsplit(as.vector(dfSel[k, "Term"]), "\\(GO"))[2]
# part2 <- paste0("\\(GO", part2)
#
# if (nchar(part1) > 40 ){
# dfSel[k, "Term"] <- paste0(part1, " \\n", part2)
# } else {
# dfSel[k, "Term"] <- paste0(part1, " ", part2)
# }
# }
# }
#
#
# #dfSel$Term <- gsub("\\(GO", "\\\n\\(GO", dfSel$Term)
#
# dfSel$Term <- factor(dfSel$Term, levels = unique(dfSel$Term))
#
# plotList[[paste0("ENR_", j)]] <- ggplot(
# data=dfSel, aes(x= Term, y=log10FDR, fill=Category, order=log10FDR)
# ) + geom_bar(stat="identity", colour="black"
# ) + coord_flip() +scale_fill_manual(values=c("yellow", "blue")) + theme(
# axis.text.y = element_text(size=8),
# axis.text.x = element_text(size=8),
# axis.title.y = element_text(size=8),
# axis.title.x = element_text(size=8),
# axis.line = element_line(colour = "black"),
# panel.border = element_rect(colour = "black", fill=NA, size=1),
# plot.title = element_text(hjust = 0.5, size = 12)
# ) + labs(title = paste0("Cluster ", clusterIDs[j]," enriched genes") ,y = "-log10(FDR)", x = ""
# ) + geom_hline(yintercept = c(-log10(0.05), log10(0.05)), color = "red", size=0.5, lty=2
# ) + geom_hline(yintercept = 0, color = "black", size=0.5
# )
# cat(" \n")
#
#
#
# ## Save to file ##
# FNbase <- paste0("Cluster_", clusterIDs[j],".enriched.genes", VersionPdfExt)
# FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
# FNrel <- paste0("report_figures/", FNbase)
#
#
# pdf(FN)
# print(plotList[[paste0("ENR_", j)]])
# dev.off()
#
# ## Create R markdown chunk ##
# link <- paste0(
# "https://",urlString,"/",
# Obio@parameterList$project_id,
# "/category-view?category_type=GO-BP"
# )
#
#
# figLegend <- paste0(
# "**Figure ",
# figureCount,
# "**: GO-BP category enrichment analysis for genes that are <font color = \\'yellow\\'>higher</font> and <font color = \\'blue\\'>lower</font> expressed in Cluster ",
# clusterIDs[j],
# " as compared to all other clusters. Download a pdf of this figure <a href='",FNrel,"' target='_blank'>here</a>. To view these gene sets in the context of your data, go to [CategoryView > GO-BP](",link,") and find the above categories using the search box."
# )
# figureCount <- figureCount + 1
#
# NewChnk <- paste0(
# "#### ", clusterIDs[j],
# "\n```{r enrchr_cluster_",
# clusterIDs[j],", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",
# figLegend,"'}\n",
# "\n",
# "\n print(plotList[['",paste0("ENR_", j),"']])",
# "\n cat( '\n')",
# "\n\n\n```\n"
# )
# }
# }
# chnkVec <- c(
# chnkVec,
# NewChnk
# )
# }###############################################################################
## Do category enrichment on clusters ##
#cat(paste(knit(text = chnkVec, quiet = T), collapse = '\n'))
## Done doing enrichment on clusters ##
###############################################################################Gene Set Heatmaps
## Two options: full heatmap and averaged heatmaps
## https://satijalab.org/seurat/v3.0/interaction_vignette.html
###############################################################################
## Add percentage expressed genes ##
# DefaultAssay(OsC) <- "RNA"
# my_genes <- rownames(x = OsC@assays$RNA)
#
# exp <- FetchData(OsC, my_genes)
#
# ExprMatrix <- round(as.matrix(colMeans(exp > 0)) *100,1)
# colnames(ExprMatrix)[1] <- "count_cut_off"
# dfExprMatrix <- data.frame(ExprMatrix)
# dfExprMatrix[["gene"]] <- row.names(dfExprMatrix)
#
# Obio@dataTableList[["dfPercCellsExpr"]] <- dfExprMatrix
#
# hmRelevantGenes <- as.vector(unique(dfExprMatrix[dfExprMatrix$count_cut_off > Obio@parameterList$singleCellPercExpressedMinCutOff, "gene"]))
#
#
#
# ## Done adding percentage expressed ##
# ###############################################################################
#
#
# ###############################################################################
# ## Make plot according to reference categories ##
# allGenes <- rownames(x = OsC@assays$RNA)
# OsC <- ScaleData(OsC, verbose = FALSE, features=allGenes)
#
# DefaultAssay(OsC) <- "RNA"
#
#
# ## Add heatmap identities to meta.data ##
#
OsC@meta.data[["hmIdent"]] <- paste0(
OsC@meta.data[,Obio@parameterList$singleCellClusterString],
"_",
substr(OsC@meta.data$sampleID,1,10)
)
#
if (length(unique(OsC@meta.data$hmIdent)) > 25){
OsC@meta.data[["hmIdent"]] <- OsC@meta.data[,Obio@parameterList$singleCellClusterString]
}
#
# Idents(OsC) <- "hmIdent"
#
# ## Done adding heatmap identities to meta.data ##
#
# printPdf <- TRUE
# referenceList <- Obio@dataTableList$referenceList
# for (i in 1:length(referenceList)){
# HMname <- names(referenceList)[i]
# geneVec <- unique(referenceList[[i]][referenceList[[i]] %in% rownames(x = OsC@assays$RNA)])
#
# if (length(geneVec) > 50){
# geneVec <- geneVec[geneVec %in% hmRelevantGenes]
# }
#
#
# ## Do Heatmap ##
# if (length(geneVec) < 1500 & length(geneVec) > 2){
# Idents(OsC) <- "hmIdent"
# HMname <- names(referenceList)[i]
# cat("\n")
# cat(paste0("**Heatmap ", HMname,"**"))
# cat("\n")
# cat("\n")
# ## Cluster genes ##
# HMgenes <- referenceList[[i]]
# #dfCluster <- OsC@assays$integrated
# Mexpr <- GetAssayData(object = OsC, assay.type = "integrated", slot = "scale.data")
# HMgenesSel <- HMgenes[HMgenes %in% row.names(Mexpr)]
#
# if (length(HMgenesSel) > 2){
# Mexpr <- Mexpr[HMgenesSel,]
#
# pdf(paste0("temp", VersionPdfExt))
# hmRes <- make.hm(
# m.df1 = Mexpr,
# filename = "",
# k.number = 1,
# n.colors = 1000,
# hclust.method = "complete",
# dist.method = "euclidean",
# main = "",
# Colv = TRUE,
# showRowNames = TRUE,
# showColNames = F,
# plotSeparationLines = FALSE
# )
# dev.off()
#
# orderedGenes <- as.vector(unique(row.names(hmRes$sorted)))
#
# if (length(unique(OsC@meta.data$hmIdent)) > 10){
# removeKey <- TRUE
# } else {
# removeKey <- FALSE
# }
#
# if (length(orderedGenes) <= 50){
# label = TRUE
# } else {
# label = FALSE
# }
#
# p1 <- DoHeatmap(
# object = OsC,
# features = orderedGenes,
# #group.by = "hmIdent",
# draw.lines =T,
# label = label,
# group.bar = TRUE,
# slot = "scale.data",
# lines.width = 2 #With of separation lines in 'cells'
# #slim.col.label = TRUE,
# #remove.key = removeKey
# ) + theme(legend.position = "none")
#
# print(p1)
#
# ## Save to file ##
# FNbase <- paste0("HM", HMname, VersionPdfExt)
# FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
# FNrel <- paste0("report_figures/", FNbase)
#
# pdf(FN)
# print(p1)
# dev.off()
#
# cat("\n")
# cat(paste0('Download a pdf of this figure [here](', FNrel, '). '))
# cat("\n")
# cat("\n")
#
# }
# }
# }
## Done making plots according to gene categories ##
###############################################################################Example Feature Plots
DefaultAssay(OsC) <- "RNA"
cat(paste0(
'Feature plots for any gene in this experiment can be viewed <a href="https://biologic.crick.ac.uk/mdata/',Obio@parameterList$project_id,'/html/FeatureView.html" target="_blank">here</a>.'
))Feature plots for any gene in this experiment can be viewed here.
plotGenes <- c(
Obio@dataTableList$referenceList$integrated_top30var[1:10]
)
plotParts <- ceiling(length(plotGenes)/2)
chnkVec <- as.vector(NULL, mode = "character")
plotListF <- list()
for (i in 1:plotParts){
tag1 <- paste0("Featureplot_",i)
featureGenes <- c(plotGenes[((2*i)-1)], plotGenes[((2*i))])
plotListF[[tag1]] <- FeaturePlot(
OsC,
features = featureGenes,
#split.by = "orig.ident",
reduction = Obio@parameterList$primReduction
)
## Save to file ##
FNbase <- paste0("Featureplot.", plotGenes[((2*i)-1)], ".", plotGenes[((2*i))],".", VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotListF[[tag1]])
dev.off()
linkFeatureView <- paste0('<a href="https://',urlString,'/mdata/',Obio@parameterList$project_id,'/html/FeatureView.html" target="_blank">FeatureView</a>')
figLegend <- paste0(
"**Figure ",
figureCount,
":** Gene expression plot for genes ",
plotGenes[((2*i)-1)],
" and ",
plotGenes[((2*i))],".",
" Results for any other gene may be plotted in ",linkFeatureView,"."
)
NewChnk <- paste0(
"#### Featureplot ", plotGenes[((2*i)-1)], " and ",plotGenes[((2*i))],
"\n```{r FeaturePlot_", i,
", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",
figLegend,"'}\n",
"\n",
"\n print(plotListF[['",tag1,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
chnkVec <- c(
chnkVec,
NewChnk
)
}cat(paste(knit(text = chnkVec, quiet = T), collapse = '\n'))Featureplot nts and alk
Figure 48: Gene expression plot for genes nts and alk. Results for any other gene may be plotted in FeatureView.
Featureplot insl5a and nrxn3a
Figure 48: Gene expression plot for genes insl5a and nrxn3a. Results for any other gene may be plotted in FeatureView.
Featureplot chgb and slc5a7a
Figure 48: Gene expression plot for genes chgb and slc5a7a. Results for any other gene may be plotted in FeatureView.
Featureplot pappaa and hs3st2
Figure 48: Gene expression plot for genes pappaa and hs3st2. Results for any other gene may be plotted in FeatureView.
Featureplot cftr and pnoca
Figure 48: Gene expression plot for genes cftr and pnoca. Results for any other gene may be plotted in FeatureView.
PCA Gene and Cell Plots
figureCount <- figureCount + 1
# An alternative heuristic method generates an 'Elbow plot': a ranking of principle components based on the percentage of variance explained by each one (`ElbowPlot` function). In this example, we can observe an 'elbow' around PC9-10, suggesting that the majority of true signal is captured in the first 10 PCs.
ElbowPlot(object = OsC) + theme(
axis.text.y = element_text(size=8),
axis.text.x = element_text(size=8),
axis.title.y = element_text(size=8),
axis.title.x = element_text(size=8),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title = element_text(hjust = 0.5, size = 12)
) + ggtitle(paste0("Variance per PCA Dimension"))
Figure, 49: Elbowplot Variance explained per PCA dimension
## Plot variance per PCA dimension #### Add PCA coordinates ##
dfTemp <- data.frame(OsC@reductions$pca@cell.embeddings)[, 1:20]
OsC <- addDf2seuratMetaData(
obj = OsC,
dfAdd = dfTemp
)
## Add UMAP coordinates to Metadata ##
dfAdd <- data.frame(OsC@reductions$umap@cell.embeddings)
OsC <- addDf2seuratMetaData(
obj = OsC,
dfAdd = dfAdd
)
## Add tSNE coordinates to Metadata ##
dfAdd <- data.frame(OsC@reductions$tsne@cell.embeddings)
OsC <- addDf2seuratMetaData(
obj = OsC,
dfAdd = dfAdd
)
Obio@dataTableList[["meta.data"]] <- OsC@meta.data
xVec <- c("PC_1","PC_3","PC_5","PC_7","PC_9","PC_11","PC_13","PC_15","PC_17","PC_19")
yVec <- c("PC_2","PC_4","PC_6","PC_8","PC_10","PC_12","PC_14","PC_16","PC_18","PC_20")
pcVec <- c("PC_1","PC_2","PC_3","PC_4","PC_5","PC_6","PC_7","PC_8","PC_9","PC_10")
chnkVec <- as.vector(NULL, mode = "character")
plotListCell <- list()
plotListGene <- list()
###############################################################################
## Collect top-enriched genes ##
EnrichedGenesList <- list()
## Done ##
###############################################################################
for (i in 1:length(xVec)){
dfDat <- Obio@dataTableList$meta.data
dfSel <- dfDat
selXY <- c(xVec[i], yVec[i])
colCol <- Obio@parameterList$singleCellClusterString
tag <- paste0(xVec[i], "and", yVec[i])
tag <- gsub("_", "", tag)
## Make Cell level PCA
plotListCell[[tag]] <- ggplot(data=dfDat, aes_string(selXY[1] , selXY[2], col=colCol, shape="sampleID")
) + geom_vline(xintercept = 0, color = "grey", size=0.5
) + geom_hline(yintercept = 0, color = "grey", size=0.5
) + geom_point()+ ggtitle(paste0("PCA - Cell Level")
) + theme(
axis.text.y = element_text(size=8),
axis.text.x = element_text(size=8),
axis.title.y = element_text(size=8),
axis.title.x = element_text(size=8),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title = element_text(hjust = 0.5, size = 12)
)
## Save to file ##
FNbase <- paste0("PCA.cell.level.", xVec[i],".", yVec[i], ".", VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotListCell[[tag]])
dev.off()
link <- paste0('<a href="https://',urlString,'/',Obio@parameterList$project_id,'/pca?x_axis=',gsub('_', '', xVec[i]),'&y_axis=',gsub('_', '', yVec[i]),'" target = "_blank">here</a>')
figCap <- paste0(
'**Figure, ' ,figureCount,'A:** Cell-level PCA plot for dimensions ', xVec[i], ' and ', yVec[i],'. ',
'Download a pdf of this figure <a href="',FNrel,'" target="_blank">here</a>. ',
'An interactive version of this figure can be found ', link, '. '
)
NewChnk <- paste0(
"#### PCA Cell Level ", xVec[i], " and ",yVec[i],
"\n```{r PCAcells_", i,
", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",
figCap,"'}\n",
"\n",
"\n print(plotListCell[['",tag,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
chnkVec <- c(
chnkVec,
NewChnk
)
## Done with cells ##
###########################################################################
###########################################################################
## Start with genes ##
dfPCADat <- data.frame(Loadings(OsC, reduction = "pca"))
dfPCADat[["gene"]] <- row.names(dfPCADat)
dfPCADat <- gather(
dfPCADat,
condition,
measurement, 1:(ncol(dfPCADat)-1),
factor_key=TRUE
)
Obio@dataTableList[["dfPCAloadings"]] <- dfPCADat
## Make Gene Level PCA ##
dfPCADat <- data.frame(Loadings(OsC, reduction = "pca"))
dfPCADat[["gene"]] <- row.names(dfPCADat)
dfPCADat <- gather(
dfPCADat,
condition,
measurement, 1:(ncol(dfPCADat)-1),
factor_key=TRUE
)
dfLoad <- dfPCADat
Obio@dataTableList$dfPCAloadings <- dfLoad
## Plot ##
selXY <- c(xVec[i], yVec[i])
dfSel <- filter(dfLoad, condition %in% selXY)
dfSel <- dfSel %>% spread(key=condition, value=measurement)
row.names(dfSel) <- dfSel$gene
dfSel[["highlight"]] <- ""
dfSel <- dfSel[order(dfSel[,selXY[1]], decreasing = FALSE), ]
dfSel[1:15, "highlight"] <- "+"
EnrichedGenesList[[paste0(selXY[1], "_neg")]]<- as.vector(dfSel$gene[1:15])
dfSel <- dfSel[order(dfSel[,selXY[1]], decreasing = TRUE), ]
dfSel[1:15, "highlight"] <- "+"
EnrichedGenesList[[paste0(selXY[1], "_pos")]]<- as.vector(dfSel$gene[1:15])
dfSel <- dfSel[order(dfSel[,selXY[2]], decreasing = FALSE), ]
dfSel[1:15, "highlight"] <- "+"
EnrichedGenesList[[paste0(selXY[2], "_neg")]]<- as.vector(dfSel$gene[1:15])
dfSel <- dfSel[order(dfSel[,selXY[2]], decreasing = TRUE), ]
dfSel[1:15, "highlight"] <- "+"
EnrichedGenesList[[paste0(selXY[2], "_pos")]]<- as.vector(dfSel$gene[1:15])
plotListGene[[tag]] <- ggplot(data=dfSel, aes_string(x=selXY[1],y=selXY[2], col="highlight")
) + geom_vline(xintercept = 0, color = "grey", size=0.5
) + geom_hline(yintercept = 0, color = "grey", size=0.5) + geom_point() + scale_color_manual(values=c("black", "red")) + ggtitle(paste0("PCA - Cell Level")
) + theme(
axis.text.y = element_text(size=8),
axis.text.x = element_text(size=8),
axis.title.y = element_text(size=8),
axis.title.x = element_text(size=8),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title = element_text(hjust = 0.5, size = 12)
)
points <- as.vector(unique(dfSel[dfSel$highlight=="+", "gene"]))
plotListGene[[tag]] <- LabelPoints(plot = plotListGene[[tag]], points =points, repel = TRUE, xnudge = 0, ynudge = 0)
## Save to file ##
FNbase <- paste0("PCA.cell.level.", xVec[i],".", yVec[i],".", VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotListGene[[tag]])
dev.off()
dim1 <- gsub("PC_", "", xVec[i])
dim2 <- gsub("PC_", "", yVec[i])
link <- paste0(
'<a href="https://',urlString,'/',
Obio@parameterList$project_id,
'/scatterplot?x_axis=add_counts_PCA_Dim_',
dim1,
'_Loadings&y_axis=add_counts_PCA_Dim_',
dim2,
'_Loadings&highlight_gene=&cat_id=ag_lab_categories__10',
'" target="_blank">here</a>'
)
figCap <- paste0(
"**Figure, " ,figureCount,"B:**Gene-level PCA plot for dimensions ", xVec[i], " and ", yVec[i], ". An interactive version of this figure can be found ", link, ". "
)
NewChnk <- paste0(
"\n```{r PCA_gene_level_", i,
", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",
figCap,"'}\n",
"\n",
"\n print(plotListGene[['",tag,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
chnkVec <- c(
chnkVec,
NewChnk
)
## Done with genes ##
###########################################################################
figureCount <- figureCount + 1
}cat(paste(knit(text = chnkVec, quiet = T), collapse = '\n'))PCA Cell Level PC_1 and PC_2
Figure, 49A: Cell-level PCA plot for dimensions PC_1 and PC_2. Download a pdf of this figure here. An interactive version of this figure can be found here.
Figure, 49B:Gene-level PCA plot for dimensions PC_1 and PC_2. An interactive version of this figure can be found here.
PCA Cell Level PC_3 and PC_4
Figure, 50A: Cell-level PCA plot for dimensions PC_3 and PC_4. Download a pdf of this figure here. An interactive version of this figure can be found here.
Figure, 50B:Gene-level PCA plot for dimensions PC_3 and PC_4. An interactive version of this figure can be found here.
PCA Cell Level PC_5 and PC_6
Figure, 51A: Cell-level PCA plot for dimensions PC_5 and PC_6. Download a pdf of this figure here. An interactive version of this figure can be found here.
Figure, 51B:Gene-level PCA plot for dimensions PC_5 and PC_6. An interactive version of this figure can be found here.
PCA Cell Level PC_7 and PC_8
Figure, 52A: Cell-level PCA plot for dimensions PC_7 and PC_8. Download a pdf of this figure here. An interactive version of this figure can be found here.
Figure, 52B:Gene-level PCA plot for dimensions PC_7 and PC_8. An interactive version of this figure can be found here.
PCA Cell Level PC_9 and PC_10
Figure, 53A: Cell-level PCA plot for dimensions PC_9 and PC_10. Download a pdf of this figure here. An interactive version of this figure can be found here.
Figure, 53B:Gene-level PCA plot for dimensions PC_9 and PC_10. An interactive version of this figure can be found here.
PCA Cell Level PC_11 and PC_12
Figure, 54A: Cell-level PCA plot for dimensions PC_11 and PC_12. Download a pdf of this figure here. An interactive version of this figure can be found here.
Figure, 54B:Gene-level PCA plot for dimensions PC_11 and PC_12. An interactive version of this figure can be found here.
PCA Cell Level PC_13 and PC_14
Figure, 55A: Cell-level PCA plot for dimensions PC_13 and PC_14. Download a pdf of this figure here. An interactive version of this figure can be found here.
Figure, 55B:Gene-level PCA plot for dimensions PC_13 and PC_14. An interactive version of this figure can be found here.
PCA Cell Level PC_15 and PC_16
Figure, 56A: Cell-level PCA plot for dimensions PC_15 and PC_16. Download a pdf of this figure here. An interactive version of this figure can be found here.
Figure, 56B:Gene-level PCA plot for dimensions PC_15 and PC_16. An interactive version of this figure can be found here.
PCA Cell Level PC_17 and PC_18
Figure, 57A: Cell-level PCA plot for dimensions PC_17 and PC_18. Download a pdf of this figure here. An interactive version of this figure can be found here.
Figure, 57B:Gene-level PCA plot for dimensions PC_17 and PC_18. An interactive version of this figure can be found here.
PCA Cell Level PC_19 and PC_20
Figure, 58A: Cell-level PCA plot for dimensions PC_19 and PC_20. Download a pdf of this figure here. An interactive version of this figure can be found here.
Figure, 58B:Gene-level PCA plot for dimensions PC_19 and PC_20. An interactive version of this figure can be found here.
Characterize PCA Dimensions
selVec <- c(Obio@parameterList$singleCellClusterString, names(OsC@meta.data)[grep("^PC", names(OsC@meta.data))])
dfPCdat <- OsC@meta.data[, selVec]
ymin <- 1.1*min(dfPCdat$PC_1)
ymax <- 1.1*max(dfPCdat$PC_1)
clusterVec <- sort(unique(OsC@meta.data[,Obio@parameterList$singleCellClusterString]))
library(scales)
clusterColVec <- hue_pal()(length(clusterVec))
i=1
chnkVec <- as.vector(NULL, mode="character")
plotList <- list()
for (i in 1:length(clusterVec)){
dfTemp <- dfPCdat[dfPCdat[,Obio@parameterList$singleCellClusterString] == clusterVec[i],]
dfTemp[,Obio@parameterList$singleCellClusterString] <- NULL
library(tidyr)
dfTemp <- gather(dfTemp, PC)
orderVec <- sort(as.numeric(gsub("PC_", "",unique(dfTemp$PC))))
orderVec <- paste0("PC_", orderVec)
dfTemp$PC <- factor(dfTemp$PC, levels = orderVec)
Ncolumns <- length(unique(dfTemp$PC))
a <- paste0("Cluster_", clusterVec[i])
tag <- paste0("PCA_Distributions_", a)
plotList[[tag]] <-ggplot(
dfTemp,
aes(x=PC, y=value, fill = PC)
) + geom_hline(yintercept = 0, color = "black", size=0.5
) + geom_jitter(width=0.1,alpha=0.2
) + geom_boxplot(
) + theme(
legend.position = "none",
axis.text.y = element_text(size=8),
axis.text.x = element_text(size=8, angle = 45,vjust = 1, hjust=1),
axis.title.y = element_text(size=8),
axis.title.x = element_text(size=8),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title = element_text(hjust = 0.5, size = 12)
) + ggtitle(paste0("PCA Distribution: ", a)
) + ylim(ymin, ymax) + scale_fill_manual(values=rep(clusterColVec[i], Ncolumns))
## Save to file ##
FNbase <- paste0(tag, VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotList[[tag]])
dev.off()
figCap <- paste0(
'**Figure, ' ,
figureCount,
':** This plot may help you to identify PCA dimensions, in which marker genes for cluster ',
clusterVec[i],
' become evident. Download a pdf of this figure <a href="',FNrel,'" target="_blank">here</a>. '
)
figureCount <- figureCount + 1
NewChnk <- paste0(
"#### ", tag,
"\n```{r ", tag,
", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",
figCap,"'}\n",
"\n",
"\n print(plotList[['",tag,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
chnkVec <- c(
chnkVec,
NewChnk
)
}cat(paste(knit(text = chnkVec, quiet = T), collapse = '\n'))PCA_Distributions_Cluster_0
Figure, 59: This plot may help you to identify PCA dimensions, in which marker genes for cluster 0 become evident. Download a pdf of this figure here.
PCA_Distributions_Cluster_1
Figure, 60: This plot may help you to identify PCA dimensions, in which marker genes for cluster 1 become evident. Download a pdf of this figure here.
PCA_Distributions_Cluster_2
Figure, 61: This plot may help you to identify PCA dimensions, in which marker genes for cluster 2 become evident. Download a pdf of this figure here.
PCA_Distributions_Cluster_3
Figure, 62: This plot may help you to identify PCA dimensions, in which marker genes for cluster 3 become evident. Download a pdf of this figure here.
PCA_Distributions_Cluster_4
Figure, 63: This plot may help you to identify PCA dimensions, in which marker genes for cluster 4 become evident. Download a pdf of this figure here.
PCA_Distributions_Cluster_5
Figure, 64: This plot may help you to identify PCA dimensions, in which marker genes for cluster 5 become evident. Download a pdf of this figure here.
PCA_Distributions_Cluster_6
Figure, 65: This plot may help you to identify PCA dimensions, in which marker genes for cluster 6 become evident. Download a pdf of this figure here.
PCA_Distributions_Cluster_7
Figure, 66: This plot may help you to identify PCA dimensions, in which marker genes for cluster 7 become evident. Download a pdf of this figure here.
PCA_Distributions_Cluster_8
Figure, 67: This plot may help you to identify PCA dimensions, in which marker genes for cluster 8 become evident. Download a pdf of this figure here.
PCA_Distributions_Cluster_9
Figure, 68: This plot may help you to identify PCA dimensions, in which marker genes for cluster 9 become evident. Download a pdf of this figure here.
PCA_Distributions_Cluster_10
Figure, 69: This plot may help you to identify PCA dimensions, in which marker genes for cluster 10 become evident. Download a pdf of this figure here.
PCA_Distributions_Cluster_11
Figure, 70: This plot may help you to identify PCA dimensions, in which marker genes for cluster 11 become evident. Download a pdf of this figure here.
PCA_Distributions_Cluster_12
Figure, 71: This plot may help you to identify PCA dimensions, in which marker genes for cluster 12 become evident. Download a pdf of this figure here.
Characterize PCA Dimensions
library(knitr)
library(ggplot2)
#save.image("temp.RData")
PCAdimensions <- paste0("PC_", 1:20)
plotList <- list()
chnkVec <- as.vector(NULL, mode = "character")
for (j in 1:length(PCAdimensions)){
posTestGeneSet <- as.vector(
unique(
EnrichedGenesList[[paste0(PCAdimensions[j], "_pos")]]
)
)
negTestGeneSet <- as.vector(
unique(
EnrichedGenesList[[paste0(PCAdimensions[j], "_neg")]]
)
)
## Get background gene set ##
#backgroundGeneVec <- row.names(OsC[["RNA"]]@counts)
if ((length(posTestGeneSet) >= 3) |(length(negTestGeneSet) >= 3)){
library(enrichR)
topMaxCat <- 10
dbs <- listEnrichrDbs()
dbs <- c("GO_Biological_Process_2017")
PosEnriched <- enrichr(posTestGeneSet, dbs)
for (i in 1:length(dbs)){
dfTemp <- PosEnriched[[dbs[i]]]
if (i ==1){
dfPosEnriched <- dfTemp
} else {
dfPosEnriched <- rbind(
dfPosEnriched,
dfTemp
)
}
}
dfPosEnriched[["log10FDR"]] <- -1*log10(dfPosEnriched$Adjusted.P.value)
dfPosEnriched <- dfPosEnriched[order(-dfPosEnriched$log10FDR),]
dfPosEnriched <- na.omit(dfPosEnriched)
## Negative Side ##
NegEnriched <- enrichr(negTestGeneSet, dbs)
for (i in 1:length(dbs)){
dfTemp <- NegEnriched[[dbs[i]]]
if (i ==1){
dfNegEnriched <- dfTemp
} else {
dfNegEnriched <- rbind(
dfNegEnriched,
dfTemp
)
}
}
dfNegEnriched[["log10FDR"]] <- -1*log10(dfNegEnriched$Adjusted.P.value)
dfNegEnriched <- dfNegEnriched[order(-dfPosEnriched$log10FDR),]
dfNegEnriched <- na.omit(dfNegEnriched)
dfNegSel <- dfNegEnriched
if (nrow(dfNegSel) > topMaxCat){
dfNegSel <- dfNegSel[1:topMaxCat,]
}
dfPosSel <- dfPosEnriched
if (nrow(dfPosSel) > topMaxCat){
dfPosSel <- dfPosSel[1:topMaxCat,]
}
if ((nrow(dfNegEnriched) > 0) | (nrow(dfPosEnriched) > 0)){
dfNegSel$log10FDR <- -1* dfNegSel$log10FDR
dfSel <- rbind(
dfNegSel,
dfPosSel
)
dfSel <- na.omit(dfSel)
dfSel <- dfSel[order(dfSel$log10FDR),]
dfSel$log10FDR <- round(dfSel$log10FDR, 2)
dfSel[["Category"]] <- ""
dfSel[dfSel$log10FDR >= 0, "Category"] <- "Enr."
dfSel[dfSel$log10FDR < 0, "Category"] <- "Depl."
for (k in 1:nrow(dfSel)){
if (nchar(dfSel[k, "Term"]) > 50 & length(grep("\\(GO", as.vector(dfSel[k, "Term"]))) > 0){
part1 <- unlist(strsplit(as.vector(dfSel[k, "Term"]), "\\(GO"))[1]
part1 <- substr(part1, 1, 45)
part2 <- unlist(strsplit(as.vector(dfSel[k, "Term"]), "\\(GO"))[2]
part2 <- paste0("\\(GO", part2)
if (nchar(part1) > 40 ){
dfSel[k, "Term"] <- paste0(part1, " \\n", part2)
} else {
dfSel[k, "Term"] <- paste0(part1, " ", part2)
}
}
}
#dfSel$Term <- gsub("\\(GO", "\\\n\\(GO", dfSel$Term)
dfSel$Term <- factor(dfSel$Term, levels = unique(dfSel$Term))
plotList[[paste0("PCA_ENR_", j)]] <- ggplot(
data=dfSel, aes(x= Term, y=log10FDR, fill=Category, order=log10FDR)
) + geom_bar(stat="identity", colour="black"
) + coord_flip() +scale_fill_manual(values=c("yellow", "blue")) + theme(
axis.text.y = element_text(size=8),
axis.text.x = element_text(size=8),
axis.title.y = element_text(size=8),
axis.title.x = element_text(size=8),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title = element_text(hjust = 0.5, size = 12)
) + labs(title = paste0("Cluster ", PCAdimensions[j]," enriched genes") ,y = "-log10(FDR)", x = ""
) + geom_hline(yintercept = c(-log10(0.05), log10(0.05)), color = "red", size=0.5, lty=2
) + geom_hline(yintercept = 0, color = "black", size=0.5
)
cat(" \n")
## Save to file ##
FNbase <- paste0("PCA_Cluster_", PCAdimensions[j],".enriched.genes", VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotList[[paste0("PCA_ENR_", j)]])
dev.off()
link <- paste0(
'<a href="https://', urlString, '/',
Obio@parameterList$project_id,
'/category-view?category_type=GO-BP" target="_blank">CategoryView > GO-BP</a>'
)
## Create R markdown chunk ##
figLegend <- paste0(
'**Figure ',
figureCount,
'**: GO-BP category enrichment analysis for the 15 genes that have <font color = "yellow"> the most positive </font> and <font color = "blue">the most negative</font> PCA loading values in dimension ',
PCAdimensions[j],
' associated with them. Download a pdf of this figure <a href="',FNrel,'" target="_blank">here</a>. To view these gene sets in the context of your data, go to ',link,' and find these categories using the search box.'
)
figureCount <- figureCount + 1
NewChnk <- paste0(
"#### ", PCAdimensions[j],
"\n```{r enrichr_",
j,", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",
figLegend,"'}\n",
"\n",
"\n print(plotList[['",paste0("PCA_ENR_", j),"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
}
}
chnkVec <- c(
chnkVec,
NewChnk
)
}###############################################################################
## Do category enrichment on clusters ##
cat(paste(knit(text = chnkVec, quiet = T), collapse = '\n'))PC_1
Figure 72: GO-BP category enrichment analysis for the 15 genes that have the most positive and the most negative PCA loading values in dimension PC_1 associated with them. Download a pdf of this figure here. To view these gene sets in the context of your data, go to CategoryView > GO-BP and find these categories using the search box.
PC_2
Figure 73: GO-BP category enrichment analysis for the 15 genes that have the most positive and the most negative PCA loading values in dimension PC_2 associated with them. Download a pdf of this figure here. To view these gene sets in the context of your data, go to CategoryView > GO-BP and find these categories using the search box.
PC_3
Figure 74: GO-BP category enrichment analysis for the 15 genes that have the most positive and the most negative PCA loading values in dimension PC_3 associated with them. Download a pdf of this figure here. To view these gene sets in the context of your data, go to CategoryView > GO-BP and find these categories using the search box.
PC_4
Figure 75: GO-BP category enrichment analysis for the 15 genes that have the most positive and the most negative PCA loading values in dimension PC_4 associated with them. Download a pdf of this figure here. To view these gene sets in the context of your data, go to CategoryView > GO-BP and find these categories using the search box.
PC_5
Figure 76: GO-BP category enrichment analysis for the 15 genes that have the most positive and the most negative PCA loading values in dimension PC_5 associated with them. Download a pdf of this figure here. To view these gene sets in the context of your data, go to CategoryView > GO-BP and find these categories using the search box.
PC_6
Figure 77: GO-BP category enrichment analysis for the 15 genes that have the most positive and the most negative PCA loading values in dimension PC_6 associated with them. Download a pdf of this figure here. To view these gene sets in the context of your data, go to CategoryView > GO-BP and find these categories using the search box.
PC_7
Figure 78: GO-BP category enrichment analysis for the 15 genes that have the most positive and the most negative PCA loading values in dimension PC_7 associated with them. Download a pdf of this figure here. To view these gene sets in the context of your data, go to CategoryView > GO-BP and find these categories using the search box.
PC_8
Figure 79: GO-BP category enrichment analysis for the 15 genes that have the most positive and the most negative PCA loading values in dimension PC_8 associated with them. Download a pdf of this figure here. To view these gene sets in the context of your data, go to CategoryView > GO-BP and find these categories using the search box.
PC_9
Figure 80: GO-BP category enrichment analysis for the 15 genes that have the most positive and the most negative PCA loading values in dimension PC_9 associated with them. Download a pdf of this figure here. To view these gene sets in the context of your data, go to CategoryView > GO-BP and find these categories using the search box.
PC_10
Figure 81: GO-BP category enrichment analysis for the 15 genes that have the most positive and the most negative PCA loading values in dimension PC_10 associated with them. Download a pdf of this figure here. To view these gene sets in the context of your data, go to CategoryView > GO-BP and find these categories using the search box.
PC_11
Figure 82: GO-BP category enrichment analysis for the 15 genes that have the most positive and the most negative PCA loading values in dimension PC_11 associated with them. Download a pdf of this figure here. To view these gene sets in the context of your data, go to CategoryView > GO-BP and find these categories using the search box.
PC_12
Figure 83: GO-BP category enrichment analysis for the 15 genes that have the most positive and the most negative PCA loading values in dimension PC_12 associated with them. Download a pdf of this figure here. To view these gene sets in the context of your data, go to CategoryView > GO-BP and find these categories using the search box.
PC_13
Figure 84: GO-BP category enrichment analysis for the 15 genes that have the most positive and the most negative PCA loading values in dimension PC_13 associated with them. Download a pdf of this figure here. To view these gene sets in the context of your data, go to CategoryView > GO-BP and find these categories using the search box.
PC_14
Figure 85: GO-BP category enrichment analysis for the 15 genes that have the most positive and the most negative PCA loading values in dimension PC_14 associated with them. Download a pdf of this figure here. To view these gene sets in the context of your data, go to CategoryView > GO-BP and find these categories using the search box.
PC_15
Figure 86: GO-BP category enrichment analysis for the 15 genes that have the most positive and the most negative PCA loading values in dimension PC_15 associated with them. Download a pdf of this figure here. To view these gene sets in the context of your data, go to CategoryView > GO-BP and find these categories using the search box.
PC_16
Figure 87: GO-BP category enrichment analysis for the 15 genes that have the most positive and the most negative PCA loading values in dimension PC_16 associated with them. Download a pdf of this figure here. To view these gene sets in the context of your data, go to CategoryView > GO-BP and find these categories using the search box.
PC_17
Figure 88: GO-BP category enrichment analysis for the 15 genes that have the most positive and the most negative PCA loading values in dimension PC_17 associated with them. Download a pdf of this figure here. To view these gene sets in the context of your data, go to CategoryView > GO-BP and find these categories using the search box.
PC_18
Figure 89: GO-BP category enrichment analysis for the 15 genes that have the most positive and the most negative PCA loading values in dimension PC_18 associated with them. Download a pdf of this figure here. To view these gene sets in the context of your data, go to CategoryView > GO-BP and find these categories using the search box.
PC_19
Figure 90: GO-BP category enrichment analysis for the 15 genes that have the most positive and the most negative PCA loading values in dimension PC_19 associated with them. Download a pdf of this figure here. To view these gene sets in the context of your data, go to CategoryView > GO-BP and find these categories using the search box.
PC_20
Figure 91: GO-BP category enrichment analysis for the 15 genes that have the most positive and the most negative PCA loading values in dimension PC_20 associated with them. Download a pdf of this figure here. To view these gene sets in the context of your data, go to CategoryView > GO-BP and find these categories using the search box.
## Done doing enrichment on clusters ##
###############################################################################Pseudotime
## Add default ##
pos <- grep("addDmaps", names(Obio@parameterList))
if (length(pos) == 0){
Obio@parameterList[["addDmaps"]] <- FALSE
}
if (Obio@parameterList$addDmaps){
library(knitr)
library(ggplot2)
invertPT <- TRUE
plotList <- list()
chnkVec <- as.vector(NULL, mode = "character")
library(Seurat)
library(destiny)
dfPCA <- OsC@meta.data
dfPCA <- dfPCA[,grep("PC_", names(dfPCA))]
dmPCA <- DiffusionMap(dfPCA)
#dpt <- DPT(dm, tips = 268)
#dpt <- DPT(dm)
dpt <- DPT(dmPCA)
#pseudotime <- dpt$dpt
# Plot DC1 vs DC2 and color the cells by their inferred diffusion pseudotime.
# We can accesss diffusion pseudotime via dpt$dpt.
df <- data.frame(
DC1 = eigenvectors(dmPCA)[, 1],
DC2 = eigenvectors(dmPCA)[, 2],
DC3 = eigenvectors(dmPCA)[, 3],
"DM_Pseudotime" = dpt$dpt
)
df$cellID <- row.names(dfPCA)
## For this project reverse pseudotime ##
if (invertPT){
PTmax <- max(df$DM_Pseudotime)
df$DM_Pseudotime <- -1* (df$DM_Pseudotime - PTmax)
## Invert DC1, 2, 3
df$DC1 <- -1* df$DC1
df$DC2 <- -1* df$DC2
df$DC3 <- -1* df$DC3
}
## Add to table ##
df$cellID <- row.names(dfPCA)
dfdbTable <- OsC@meta.data
dfdbTable$row_names <- NULL
dfdbTable$DC1 <- NULL
dfdbTable$DC3 <- NULL
dfdbTable$DC2 <- NULL
dfdbTable$DM_Pseudotime <- NULL
dfdbTable <- unique(dfdbTable)
###############################################################################
## Add pseudotime components ##
#dim(dfdbTable)
dfdbTable[["cellID"]] <- row.names(dfdbTable)
dfdbTable <- merge(
dfdbTable,
df,
by.x = "cellID",
by.y = "cellID"
)
#dim(dfdbTable)
## Add to Seurat object ##
OsC@meta.data <- dfdbTable
## Create Pseudotime plot ##
dfTemp <- dfdbTable
#dotsize <- 0.5
dotcolor <- "darkblue"
tag <- "PC1PC2all"
plotList[[tag]] <- ggplot(dfTemp, aes(PC_1, PC_2, color=DM_Pseudotime)
)+ geom_point(
shape = 16,
size = as.numeric(dotsize)
) + xlab("PC1") + ylab("PC2") + scale_color_gradient2(
low="#ff6600",
high=dotcolor #,
#limits=c(0,maxExpr)
) + theme(
axis.text.y = element_text(size=8),
axis.text.x = element_text(size=8),
axis.title.y = element_text(size=8),
axis.title.x = element_text(size=8),
axis.line = element_line(colour = "black"),
panel.border = element_rect(colour = "black", fill=NA, size=1),
plot.title = element_text(hjust = 0.5, size = 12),
panel.background = element_rect(fill = "lightgrey")
) + ggtitle("PC1, PC2 and DM Pseudotime"
) #+ xlim(minX, maxX) + ylim(minY, maxY)
## Save to file ##
FNbase <- paste0("Pseudotime_overview", VersionPdfExt)
FN <- paste0(Obio@parameterList$reportFigDir, FNbase)
FNrel <- paste0("report_figures/", FNbase)
pdf(FN)
print(plotList[[tag]])
dev.off()
## Create R markdown chunk ##
figLegend <- paste0(
'**Figure ',
figureCount,
'**: Figure depicting PCA components 1 and 2 with the diffusion map pseudotime highlighted in color. Download a pdf of this figure <a href="',FNrel,'" target="_blank">here</a>. '
)
figureCount <- figureCount + 1
NewChnk <- paste0(
"\n#### Pseudotime All Timepoints",
"\n```{r ", tag, ", results='asis', echo=F, eval=TRUE, warning=FALSE, fig.cap='",
figLegend,"'}\n",
"\n",
"\n print(plotList[['",tag,"']])",
"\n cat( '\n')",
"\n\n\n```\n"
)
chnkVec <- c(
chnkVec,
NewChnk
)
} # end dmap###############################################################################
## Do category enrichment on clusters ##
if (Obio@parameterList$addDmaps){
cat(paste(knit(text = chnkVec, quiet = T), collapse = '\n'))
}Pseudotime All Timepoints
Figure 92: Figure depicting PCA components 1 and 2 with the diffusion map pseudotime highlighted in color. Download a pdf of this figure here.
## Done doing enrichment on clusters ##
################################################################################## Will save Obio object here, so it can be re-used with different parameters
save(Obio,
file = paste0(
Obio@parameterList$localWorkDir,
Obio@parameterList$project_id,
".bioLOGIC.Robj"
)
)
print("Obio Object saved.")
save(OsC,
file = paste0(
Obio@parameterList$localWorkDir,
Obio@parameterList$project_id,
".Seurat.Robj"
)
)## Try to retrieve project data from db ##
library(RMySQL)
db.pwd2 <- "_asf_"
db.user2 <- "asf"
host2 <- "ms1.thecrick.org"
projectParams <- Obio@documentationParams
tryCatch({
dbDB = dbConnect(drv = RMySQL::MySQL(), user = db.user2, password = db.pwd2, host = host2, dbname = "asf");
dfProposal = dbGetQuery(dbDB, paste0("SELECT * FROM asf_proposals WHERE project_name ='",Obio@parameterList$lims.id,"'"));
dbDisconnect(dbDB)
}, error = function(x) {
message("Project Database could not be reached or has no entry in Obio@parameterList$lims.id for this analysis.")
})[1] TRUE
if (nrow(dfProposal) == 1){
if (!is.na(dfProposal[1,"ProjectAlias"]) & dfProposal[1,"ProjectAlias"] != ""){
projectParams[["title"]] = paste0(dfProposal[1,"ProjectAlias"], " - ", dfProposal[1,"project_name"])
}
if (!is.na(dfProposal[1,"project_user"]) & dfProposal[1,"project_user"] != ""){
projectParams[["subtitle"]] = paste0(dfProposal[1,"user_lab"], " Lab - ", dfProposal[1,"project_user"])
projectParams[["subtitle"]] <- gsub("^ Lab - ", "", projectParams[["subtitle"]])
}
if (!is.na(dfProposal[1,"proposal_text"]) & dfProposal[1,"proposal_text"] != ""){
projectParams[["abstract"]] = dfProposal[1,"proposal_text"]
}
}
## Escape all special characters
projectParams <- lapply(
projectParams, function(x)
#gsub("([.|()\\^{}+$*?]|\\[|\\])", "\\\\\1", x)
gsub("([.|()/\\^{}+$*?]|\\[|\\])", " ", x)
)
projectParams <- lapply(
projectParams, function(x)
#gsub("([.|()\\^{}+$*?]|\\[|\\])", "\\\\\1", x)
gsub("\\\n", " ", x)
)
#projectParams$title <- "Title"
# projectParams$abstract <- "This is the QC section."
#projectParams$subtitle <- "Abstract"Documentation
sessionInfo()## R version 3.6.0 (2019-04-26)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: CentOS Linux 7 (Core)
##
## Matrix products: default
## BLAS: /camp/apps/misc/stp/babs/manual/software/r/R-3.6.0-foss-2016b/lib64/R/lib/libRblas.so
## LAPACK: /camp/apps/misc/stp/babs/manual/software/r/R-3.6.0-foss-2016b/lib64/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_GB.utf-8 LC_NUMERIC=C
## [3] LC_TIME=en_GB.utf-8 LC_COLLATE=en_GB.utf-8
## [5] LC_MONETARY=en_GB.utf-8 LC_MESSAGES=en_GB.utf-8
## [7] LC_PAPER=en_GB.utf-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_GB.utf-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] gplots_3.0.3 RColorBrewer_1.1-2
## [3] yaml_2.2.1 openxlsx_4.1.5
## [5] glmGamPoi_1.1.13 destiny_2.14.0
## [7] enrichR_2.1 DT_0.14
## [9] AUCell_1.6.1 sleepwalk_0.3.0
## [11] ggtree_1.16.6 scales_1.1.1
## [13] ggrepel_0.8.2 RMySQL_0.10.20
## [15] DBI_1.1.0 DESeq2_1.24.0
## [17] SummarizedExperiment_1.14.1 DelayedArray_0.10.0
## [19] BiocParallel_1.18.1 matrixStats_0.56.0
## [21] Biobase_2.44.0 GenomicRanges_1.36.1
## [23] GenomeInfoDb_1.20.0 IRanges_2.18.3
## [25] S4Vectors_0.22.1 BiocGenerics_0.30.0
## [27] knitr_1.29 forcats_0.5.0
## [29] stringr_1.4.0 purrr_0.3.4
## [31] readr_1.3.1 tidyr_1.1.0
## [33] tibble_3.0.1 tidyverse_1.3.0
## [35] ggplot2_3.3.2 Seurat_3.1.5
## [37] dplyr_1.0.0
##
## loaded via a namespace (and not attached):
## [1] rappdirs_0.3.1 ggthemes_4.2.0 R.methodsS3_1.8.0
## [4] jrc_0.3.1 acepack_1.4.1 bit64_0.9-7
## [7] irlba_2.3.3 R.utils_2.9.2 data.table_1.12.8
## [10] rpart_4.1-15 RCurl_1.98-1.2 generics_0.0.2
## [13] cowplot_1.0.0 RSQLite_2.2.0 RANN_2.6.1
## [16] proxy_0.4-24 future_1.17.0 bit_1.1-15.2
## [19] xml2_1.3.2 lubridate_1.7.9 httpuv_1.5.4
## [22] assertthat_0.2.1 xfun_0.15 hms_0.5.3
## [25] evaluate_0.14 promises_1.1.1 DEoptimR_1.0-8
## [28] fansi_0.4.1 caTools_1.18.0 dbplyr_1.4.4
## [31] readxl_1.3.1 igraph_1.2.5 geneplotter_1.62.0
## [34] htmlwidgets_1.5.1 ellipsis_0.3.1 RSpectra_0.16-0
## [37] crosstalk_1.1.0.1 backports_1.1.8 annotate_1.62.0
## [40] vctrs_0.3.1 TTR_0.23-6 ROCR_1.0-11
## [43] abind_1.4-5 RcppEigen_0.3.3.7.0 withr_2.2.0
## [46] robustbase_0.93-6 checkmate_2.0.0 vcd_1.4-7
## [49] sctransform_0.2.1 treeio_1.8.2 xts_0.12-0
## [52] cluster_2.0.8 ape_5.4 segmented_1.2-0
## [55] lazyeval_0.2.2 laeken_0.5.1 crayon_1.3.4
## [58] genefilter_1.66.0 pkgconfig_2.0.3 labeling_0.3
## [61] nlme_3.1-139 nnet_7.3-12 rlang_0.4.6
## [64] globals_0.12.5 lifecycle_0.2.0 modelr_0.1.8
## [67] rsvd_1.0.3 cellranger_1.1.0 lmtest_0.9-37
## [70] graph_1.62.0 Matrix_1.2-18 carData_3.0-4
## [73] boot_1.3-22 zoo_1.8-8 reprex_0.3.0
## [76] base64enc_0.1-3 ggridges_0.5.2 png_0.1-7
## [79] viridisLite_0.3.0 rjson_0.2.20 bitops_1.0-6
## [82] R.oo_1.23.0 KernSmooth_2.23-15 DelayedMatrixStats_1.6.1
## [85] blob_1.2.1 jpeg_0.1-8.1 memoise_1.1.0
## [88] GSEABase_1.46.0 magrittr_1.5 plyr_1.8.6
## [91] ica_1.0-2 gdata_2.18.0 zlibbioc_1.30.0
## [94] compiler_3.6.0 fitdistrplus_1.1-1 cli_2.0.2
## [97] XVector_0.24.0 listenv_0.8.0 patchwork_1.0.1
## [100] pbapply_1.4-2 htmlTable_2.0.0 Formula_1.2-3
## [103] MASS_7.3-51.4 mgcv_1.8-28 tidyselect_1.1.0
## [106] stringi_1.4.6 highr_0.8 locfit_1.5-9.4
## [109] latticeExtra_0.6-29 grid_3.6.0 tools_3.6.0
## [112] future.apply_1.6.0 rio_0.5.16 rstudioapi_0.11
## [115] foreign_0.8-71 gridExtra_2.3 smoother_1.1
## [118] scatterplot3d_0.3-41 farver_2.0.3 Rtsne_0.15
## [121] digest_0.6.25 rvcheck_0.1.8 BiocManager_1.30.10
## [124] shiny_1.5.0 Rcpp_1.0.4.6 car_3.0-8
## [127] broom_0.5.6 later_1.1.0.1 RcppAnnoy_0.0.16
## [130] httr_1.4.1 AnnotationDbi_1.46.1 kernlab_0.9-29
## [133] colorspace_1.4-1 rvest_0.3.5 XML_3.99-0.3
## [136] fs_1.4.2 ranger_0.12.1 reticulate_1.16
## [139] splines_3.6.0 uwot_0.1.8 tidytree_0.3.3
## [142] sp_1.4-2 plotly_4.9.2.1 xtable_1.8-4
## [145] jsonlite_1.7.0 R6_2.4.1 Hmisc_4.4-0
## [148] pillar_1.4.4 htmltools_0.5.0 mime_0.9
## [151] glue_1.4.1 fastmap_1.0.1 VIM_6.0.0
## [154] class_7.3-15 codetools_0.2-16 utf8_1.1.4
## [157] tsne_0.1-3 lattice_0.20-38 mixtools_1.2.0
## [160] curl_4.3 leiden_0.3.3 gtools_3.8.2
## [163] zip_2.0.4 survival_3.2-3 rmarkdown_2.3
## [166] munsell_0.5.0 e1071_1.7-2 GenomeInfoDbData_1.2.1
## [169] haven_2.3.1 reshape2_1.4.4 gtable_0.3.0The Francis Crick Institute, stefan.boeing@crick.ac.uk↩︎